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ABSTRACT 


Solvent flooding is the basis of a wide range of enhanced oil 
recovery methods, and has been shown to be a possible method of 
creating initial steam injectivity in tar sands. This investigation 
presents a unique model for the dissolution of a semi-solid bitumen, 
resulting from the injection of a solvent. The solution of the 
mathematical equations, describing this phenomena is discussed. Results, 
for a series of two-dimensional problems, are presented. Numerical 
aspects are addressed. 

The mathematical treatment is based upon the rapid dispersion 
of solvent in the pores, which contain bitumen particles dissolving in 
a streaming interstitial solution. Macroscopic transport mechanisms, 
based upon the microscopic phenomena, trace the course of the diffusing 
species in a dynamically interactive solid-liquid network. The fluid 
velocity field is implicitly dependent on both solvent concentration 
and liquid saturation, while an apparent diffusion coefficient tensor 
includes the combined effects of molecular diffusion and directional 
velocity dispersion. The bitumen particles are given spherical symmetry 
such that a radial dissolution integro-differential equation accounts 
for the volume changes from solvent transfer into the initially immobile 
bitumen. 

An overall numerical approach to the nonlinear system of 
equations determines the state vector of the fluid pressure, the fluid 
phase solvent concentration and the fluid saturation. The instantaneous 
radii of the bitumen particles are determined from an analytical solution 


of the moving boundary problem associated with the transfer of solvent 


bho beomauley to ogicy obbw§ Ye pieed ade at aeitoas? snoeton 
| yo Hosgan olds 2a0q sol ad ot trois ve ih fie chosen © 
aeiregiyeovtr etal . ables cet wh. cpivisoopet nara cabs i . 
_ rsmtari d bbide sce 5 tO Anedlithan th ont. rot! ksthone seeping 2m 
ens To righ eee AT Tew o4 ‘s is pc reoe tot sda woes ; | 
ation? Oo Buak tb 7. actaaorere ety eee, aasiivenps iaaia 
fest omul panied a7 9To. peer} ighiontiee +o asi 


\ 


notetaqeth bigqad 97l4-moqy, hedged ak Shgatets, pameta 


a gaty tozetb eglors rag nesninttid nin PaeD sg fe waa ot. th 3 
Abelrodsem s1ogerert D idoyv2oge 46! rife pmlalaetuoond ga 
gavess2ibh Sit to S202 sy o3nes ‘aeowiorss le 
btult ol sows ait hrpbs 7b Eo avis onsataAt yt Pc etree a ne 
KorsteTINpoNoD Tov Loe Asoxd fOr, seadhasieesh ‘koko s fest of bigs? 

» ToRwed: . ts toe ee oukeuh> hb ney a oo bth Mh Ia Iae sae 
(aMoitsorih bits gohewyith -tetaioa tam $4. 24y1 Dortivmes sit sabud ant 


Crisis Eaotieiqe nevig ave eolattiug gegagid eff . moles it Cae 


symone 08 Sandpe inbregretPtS—onys it hotiinrsib Tether s tars dou | 


9 t Demi iasstat ost Ocirt awicnees Jueviaa wow? evgnerio amatie 


oy ae 


© mevave Tats Ihe od? oF enna feakveee Liavava BA * 
Ofer? ods .crieeenq Stett one 26 T6honv afere enr 2ankexsreb enokgaupe 7 
avouneinateni aT °.ieicormdee belt ot baie soisaronepas Jaev (oe oem 


moivuioe isstsylens 2a woxt honlersteb wae estubsreg romss te ents Yo wae ic 
tneeioe to totensrt odd 0 kee fasutooses, mad dere vrebquod given axiz ‘To. 


ae 


° & 


5 ise Pus : di 
SS J 7 - fie 


_ 


into a single particle from the streaming solution. The macroscopic 
saturation follows from a localized average particle radius. Together 
with the combined Darcy-continuity and convective-diffusion with 
adsorption equations, the coupled system is solved by a direct approach, 
using finite difference techniques. 

It is shown that three dimensionless groups control the bitumen 
leaching process. In the absence of fluid flow, the dissolution of 
bitumen is governed by the Damkéhler and solvent capacity numbers. In 
the absence of dissolution, the Peclet number governs the miscible 
displacement of an initial bitumen distribution in the fluid phase. 

This work presents computational results for a series of two- 
dimensional problems with focus on the five spot flow geometry, showing 


the importance of the model developed. 
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CHAPTER ONE 


INTRODUCTION 


Man's use of bitumen dates back to the roots of civilization, 
historically recorded as early as the legendary Tower of Babel! . 
Through the ages, man has grown accustomed to depend on petroleum. 

The Eygptians lined their graineries with bitumen and used its 
derivatives to mummify their Pharaohs, while Roman ladies applied it 

as cosmetic by the light of their petroleum fueled lamps. Present 

day civilization absorbs it insatiably while science labours ad infinitum 
to produce it. 

Mega-projects have been proposed for producing supergiant 
underground fields of tar, which have remained untouched for millennia. 
Modern technology must first overcome virtual immobility before these 
ancient fields flow. The subject of this work concerns itself with 


bitumen leaching from these massive sands as a viable recovery 


procedure with the technology at hand. 


1.1." Statement of the Problem 

Bitumen recovery from tar sands is limited by the low injectivity 
associated with the displacement of the highly Viscous tar, Miscipie, OF 
solvent flooding dilutes the oil in place and effectively increases the 
fluid mobility. 

Favourable viscosity ratios provide for the stable advancement 
of the displacement front, with the limiting solvent-oil contact surface 


propagating in a piston-like fashion. Adverse viscosity ratios initiate 
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frontal instability, with the viscous fingering of solvent by-passing 
a major portion of the oil in place. For the extreme viscosity ratios 
characteristic of tar sand, the solvent rapidly fingers through the 
medium to give an early breakthrough. Once these high mobility fingers 
have set their paths of least resistance, dissolution becomes the 
controlling mechanism for bitumen extraction. 

The mechanisms by which solvent fingers through the medium is 
beyond the scope of this work. The current problem concentrates on 
the transport phenomena after breakthrough. Earlier experiments by 
S.M. Farouq Ali have used gas injection techniques to obtain uniformly 
flooded tar sand packs with solvent at low saturations. This study 
considers the breakthrough time to be negligible with respect to the 
recovery time-scale and the injected solvent is assumed to instantaneously 
invade all interstices to give an initial continuity throughout. 

The main objectives of this work are to develop a two-dimensional 
mathematical model for the in-situ isothermal solvent leaching of 
bitumen from tar sand, the implementation of solution and the analysis of 
Simulated results. 

Macroscopic transport phenomena are used to predict the state 
vector of fluid pressure, fluid solvent concentration and fluid fraction— 
essentially describing the leaching process. Simulation provides the 
means for the numerical advancement of a mathematical flood through a 
discrete representation of the physical system, with similarity groups 
establishing the scaling criterion between reservoir and laboratory 


conditions. 
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Two mechanisms compete for leaching status — miscible displace- 
ment and dissolution kinetics. These special cases are to be investi- 
gated in order to fully understand the individual mechanisms, and their 


contributions to the overall process. 


1.2. Literature Survey 
Current literature on leaching phenomena assume small porosity 
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change in the reaction zone limited to the solution mining of 


minerals. Solution mining considers the dissolution of a porous rock 


> 


matrix by surface reaction’ » Whereas the leaching of tar sands dissolves 
the interstitial bitumen. The present work relates porosity change to 
the dissolution of spherical bitumen particles, embedded in an inert 
rock matrix, subject to solvent adsorption. 

Diffusion in spheres and linear moving boundary problems have 
been studied in great detatn >’! Ready and Cooper® modelled the 
dissolution of a single sphere in an infinite medium, accounting for 
the effects of the moving boundary and the resulting convective transport 
from the propagating interface. They assumed the interface reaction to 
be virtually instantaneous, with a constant solute concentration at the 
interface and interdiffusion as the only mechanism controlling the rate 
of interface motion. Transport within the sphere's interior was assumed 
sufficiently rapid to give a constant internal solute concentration, and 
the concentration in the surrounding solution varied with alga? 
distance only. Their main concern was to model the physical problem of 


bubble growth and dissolution. In contrast, this work assumes dispersion 


in a finite liquid phase to be sufficiently large, such that a local 
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macroscopic average represents the microscopic concentration at every 
point in the solution about a multitude of spherical tar particles. In 
particular, radial diffusion into the solid interior is included, 
accounting for the solvent stored internally by absorption. 

The solution process in polymer-solvent systems produce a thin 
transition layer”, separating solid from liquid. Rosner and Epstein?” 
investigated the effect of interface kinetics on bubble growth by 
introducing a thin concentration boundary layer to obtain closed form 
solutions. The basic concept of a thin transition layer is utilized in 
this work to relate first-order adsorption to the solvent transfer 
across the free boundary. 

Mathematically, a radial spherical moving boundary problem 
determines the instantaneous radii of dissolving particles based on an 
average local size distribution for a given number density. Material 
balance, relating absorption kinetics to the rate of solvent transfer 
across the free boundary, is described by an integro-differential 
process. The liquid pore space available for miscible displacement 
follow directly from the average radius of the spherical dissolution 
front. 

Miscible displacement is modelled by the convective-diffusion and 
Darcy-continuity isch as acim Analytical solutions are not available 
for the general miscible displacement problem due to reservoir hetero- 
geneities and the nonlinearity of the governing equations. ‘Several 


methods + = 


have been proposed for the numerical solution, but 
or ae 12 LUBA 
difficulties such as overshoot ~, numerical dispersion and grid 


orientation éffects” limit their range of applicability. 
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The Darcy-continuity equation is a second-order elliptic partial 
differential equation whereas the convective-diffusion equation behaves 
like either a second-order parabolic equation when dispersion dominates, 
or a first-order hyperbolic equation when transport dominates. 

The second-order central difference approximation developed for 
parabolic equations suffers from nonphysical oscillations, or over- 
shoot ,!7 when convection becomes significant, and a finer grid mesh 
becomes necessary to resolve concentration profiles. 1° First-order 
upstream differences for approximating the convection term eliminates 
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overshoot but introduces large numerical dispersion errors 
may be several orders of magnitude larger than the physical dispersion 
term. Implementation of these methods are advantageously simple and 
extension to multidimensional problems follow with relative ease. Their 
two-dimensional versions comprise the standard five-point difference 
formulas, serving as basis for many schemes offering an improvement in 
accuracy. 

Laumbach?’ shows that a truncation cancellation procedure using 
weighted spatial averages with central differences, gives a high-order 
method eliminating overshoot and numerical dispersion for select time- 
steps. The method is conditionally stable, and an extension to two- 
dimensional problems is outlined using the alternating direction implicit 
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ry oe uses explicit, upstream-weighted finite differences 


Larson 
in space and time with variably timed flux updating. The convected 
flux is updated with a frequency determined by the concentration 


velocity — a linear combination of a downstream-weighted convection 
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velocity with a central-differenced dispersion velocity. Numerical 
dispersion is reduced, however frontal profiles suffer from the inherent 
inaccuracy of finite-difference pressure and flow fields. 

Investigators including Settari, Price and Dupont! and Younnee 
have shown that finite elements, based on Galerkin Methods seas have 
higher accuracy than standard finite differences. 

Garder, Peaceman and Pozzi-> proposed the method of characteristics 
used for hyperbolic equations to resolve sharp fronts at high convection 
rates. A system of ordinary differential equations determine the 
position and concentration of moving points using the finite-difference 
velocities interpolated on a stationary grid. 

Ewing, Russel and Wheelen-© have considered using mixed methods 
with a modified method of characteristics. The base method utilizes 
characteristic flow directions to model convection and finite elements 
for diffusion and dispersion. They suggest that care must be exercised 
when treating Sue ee around wells and across flowing boundaries. 
The method generates some overshoot, reflecting the inability to resolve 
sharp fronts on a course mesh. 

Other schemes have been proposed using a somewhat different 
approach to the convective-diffusion problem. 

The random choice method-” is based on operator splitting with 
fractional steps treating convection and diffusion separately. The 
first fractional step solves the convection equation without diffusion 


28,29 to calculate the 


using Glimm's (piecewise sampling) method 
convective flow, and a second fractional step accounts for the parabolic 


character of the diffusion term using central differences. The method 
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resolves sharp fronts, but has first-order random errors in frontal 
position due to the sampling procedure used. 
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Schemes, including the operator compact implicit method 


2 ‘ 
233 use constraints based on the 


and single cell high order formulas,” 
highest order accuracy possible on the smallest difference star, to 
determine suitable weighting coefficients. The single cell high order 
scheme, based on a local series solution, is intrinsically multi- 
dimensional and uses a nine-point formula for the two-dimensional case. 
The operator compact implicit procedure derives a spatial operator for 
the combined convective and diffusive terms using Taylor series. 
Restricted invertibility’© of the compact operator limits this methods 
range of applicability and the extension to multidimensional problems 
involve the explicit evaluation of inverse operators for an alternating 
direction implicit procedure. 

The fundamental problem with all multidimensional numerical 
methods, and in particular five-point schemes, is the nonuniqueness of 
solution for different orientations of the grid with respect to the 
pattern flow lines. The nonuniqueness was first observed by Todd et 
See when the solution of an identical problem on two extreme grid 
orientations produced conflicting results. They proposed the two-point 
upstream mobility approximation to alleviate the grid orientation 
problem, but were only successful for mildly adverse conditions. 

Yanosik and McCracken” postulated that the grid orientation 
problem was caused by the lack of flow from a grid block to its diagonal 


neighbours, and introduced a nine-point formulation to account for this 


flow. The method alleviates grid orientation effects at moderately 
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adverse conditions, however frontal instabilities develop in both 
orientations at higher mobility ratios. 

Robertson and Woo?” proposed the use of orthogonal curvilinear 
coordinates which coincide with the streamlines and equipotentials for 
pattern floods having a unit mobility ratio, in place of conventional 
rectilinear grids, to reduce the grid orientation effect. 

Vinsome and Au>® discussed the directional dependence of the 
pressure field truncation error as the source for these orientation 
effects and suggested the use of total harmonic mobilities in place of 
standard upstream mobilities. 

Ewing et Bie s also mentioned that the numerical dispersion 
generated by truncation error is not rotationally invariant and this 
directional dependence causes grid orientation effects to occur in the 
solution. 

For the extremely adverse mobility ratios associated with the 
miscible displacement of tar, all the mentioned techniques will exhibit 
some numerical deficiency. In the combined miscible displacement and 
dissolution problem, numerical difficulties are reduced by assuming 
solvent instantaneously invades the medium to give a uniform initial 
concentration. The oil-solvent discontinuity is implicit in the 
dissolution equation and concentration gradients remain small at the 
expense of a time dependent fluid saturation. 

The solution methodology followed in this work is based on its 
simplicity of implementation. A five-point scheme using finite 
differences is the skeletal framework. A novel approach to eliminate 


numerical dispersion, reduce overshoot and alleviate grid orientation 
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effects is attempted for two-dimensional miscible floods. Extending 
Laumbach's*” work to be used in conjunction with direct solution 
algorithms, this approach transforms the third-order truncation into 
local velocity or streamline coordinates, cancelling a dominant 


longitudinal error flux term at every mesh point. 


1.3. Tensor Notation 

This work has adopted Einstein's notation convention for 
differential expressions, where summations are to be performed over 
repeated indices in order to evaluate the tensor products. Standard 
vector notation is also used liberally, with "tilde" (%) asa 


superscript to distinguish tensor quantities. 


1.4. Finite-Difference Framework 

Three approaches may be considered in the derivation of finite 
difference analogs of the differential system. The partial differential 
equations can be discretized using finite differences in which reflection 
conditions with the introduction of fictitious points outside the 
physical grid are used to satisfy the boundary conditions. Secondly 
integral finite differences are obtained by spatially discretizing the 
integral conservation laws directly, bypassing the intermediate step 
of writing down the partial differential equations. Discretization of 
the integral forms yield transmissibilities between adjacent points, 
which are readily interpreted as the conductivities at the bounding 
faces of the volume element. The use of half and quarter blocks 


eliminates the need for fictitious points and the boundary conditions 
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are satisfied by constraining the appropriate transmissibility. Where 
central-differences are concerned, both approaches invariably lead to 
the same algebraic equations. The third approach constructs finite 
difference analogs of the differential system on the basis of existing 
mesh points in the grid. To obtain boundary difference equations of 
the same order, the boundary difference star must extend deep into the 
interior to determine the equation for the specified condition at the 
peripheral point. 

The governing differential equations follow from integral 
conservation laws, and without any loss of generality, central- 
differencing techniques are applied to the integral forms. The author 
feels that integral discretization gives a better correspondence to the 
physical flow terms, and an increased compactness of notation, with 
superscripts denoting time levels and subscripts giving the spatial 


coordinates. 


ties LNeSAS Otructure 

This introductory chapter has set the groundwork for the 
following chapters. The problem has been posed, a skeletal background 
presented and a solution methodology identified. 

Chapter Two translates the physical aspects of the problem into 
mathematical language, carefully defining mechanistic dynamics with 
elaborate theoretical formulation. 

Chapter Three outlines the solution, with discrete mathematical 
algorithms aimed at high-speed digital computers. 


Chapter Four isolates the special case of miscible displacement, 
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scrutinizes simulator performance, and investigates similarity groups. 

Chapter Five studies the effect of dissolution kinetics on phase 
behaviour, combines it with miscible displacement to analyze predicted 
leaching production data. 

The final chapter draws conclusions on the basis of completed 
work, with remarks directed for further study. 

Four appendices detail rigorous mathematical derivations, 


covering analytical solutions and high-order finite difference schemes. 
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CHAPTER TWO 


THEORETICAL FORMULATION 


The mathematical derivation to the physical problem of in-situ 
isothermal solvent leaching of bitumen from tar sand is presented. The 
particular field geometry of a two-dimensional rectangular lamina with 
a two-component solvent-bitumen system is considered. A dissolution 
model for bitumen spheres is proposed. 

A series of mass balances are executed on both macroscopic and 
miscroscopic scales to derive the Darcy-continuity equation for fluid 
flow, the convective-diffusion with adsorption equation for solvent 
transport in the fluid stream, and the radial dissolution equation for 
the relative volume fractions of solid and liquid phases. Conditions 


for inlet, outlet and impermeable boundaries are established. 


2.1. The Porous Region 

Let F{Xx,} represent a stationary reference frame with the 
cartesian basis vectors {é,}. Consider a finite porous lamina, R, of 
uniform thickness, Az, whose curvature is determined by the height, h, 
as measured from a horizontal reference. The projection of R_ on the 
horizontal is a rectangle with dimension (L,W). The region R is 
bounded by a closed surface, a in which an inert rock matrix provides 
pore space for both liquid and semi-solid phases to coexist. The liquid 
phase consists of a quasi-continuous streaming solution, where constituent 
species propagate by convective and diffusive transport mechanisms. The 


semi-solid phase is composed of a discrete set of spherical solute 
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particles confined to the pores of the matrix rock. Constituent semi- 
solid species propagate by diffusion only. By definition, the streaming 
solution separates all solute particles from the rock matrix to provide 
maximum surface area for interphase mass transfer. 

Consider a finite volume element, AV, bounded by the closed 
surface, o, about a macroscopic point P(x). The instananeous porosity, 
¢, 1s defined as the volume fraction in which convection currents exist — 
the liquid phase. The time independent matrix porosity, on is the 
volume fraction spanned by the combined solid-liquid system, with the 


Saturation, 5, as the laquid: fraction. 


S = Ravi) 


Los 
®m 
Let R(x) represent the macroscopic solid subregion in the 
interior of the volume element AV(X). The solute particles have a 
uniform local size distribution such that a particle number density, 
N_, relates the average particle radius, tT)? to the macroscopic liquid 


saturation. 


2.1 (2) 


2.2. Radial Dissolution 

Concentration gradients are established at the solute-solution 
contact surface of each particle as a result of the discontinuities 
in diffusive properties. These gradients create diffusion currents on 


the molecular level, carrying bitumen from the solute particles into 
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the streaming interstitial solution. As bitumen desorbs from the 
solute, the depletion in solid mass sets the bounding spherical 
surface into motion to balance the consequent shrinkage in solid 
volume. The radial advance velocity of the solute-solution contact 


surface is defined as the velocity of dissolution, x 
ro=—2 geet) 


A fluid source instantaneously invades the medium to give an 
initial liquid saturation with uniform solvent composition throughout. 
All solid particles are initially bitumen and void of solvent. A 
critical concentration, C*, determines the solidity of the solid- 
liquid interface. By definition, the solvent concentration within the 
solid interior is less than or equal to the critical and any point 
external in the liquid phase is greater. The interfacial concentration 
buildup to critical is instantaneous and the single-valued surface 
concentration is C* for the duration of the entire dissolution 
process? 

All particles have homogeneous, isotropic diffusive properties, 
such that only normal concentration gradients exist at the contact 
surface. A semi-solid molecular diffusion coefficient, D.; determines 
the diffusion rate of solvent in the solid interior. Microscopic 
mixing in the interstitial liquid is rapid and an ambient concentration 
stabilizes outside a thin transition layer of thickness, 6*. The 
radial microscopic solvent concentration, C.; is continuous across the 


transition zone with discontinuities in slope marked by the rapid change 
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in diffusive properties across the solid-liquid contact surface as 
depicted in Figure 2.2-1. The concentration is assumed linear across 
the transition layer and the stabilized ambient approximated by the 
local macroscopic average, C. 

Diffusion is Fickian. An apparent convection flux across the 
transition layer results from the moving concentration surface. A 
solvent balance across the transition layer equates the solvent influx 
from the liquid to the flux entering the solid, 


aC , 9) siD* 
eee + C*r = ce (G2C2 cs cr, ‘ap oy 


cag 8 Pa Pp 
JOE 
with Ds as the mean radial diffusion coefficient across the transition 


zone. The adsorption coefficient, a, 


o= 2.2 (3) 


governs the rate of solvent adsorption across the transition zone. 
Conservation of solvent mass over the solid subregion equates the rate 
of solvent transfer from the liquid phase, through the transition 


layer, to the net rate of accumulation within the solid spheres. 
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Each sphere has an initial radius, To» and is subject to radial 


diffusion with a time dependent or moving boundary. 
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Figure 2.2-1 


Radial Dissolution 


U7 


pf ' 


) : : 


- 
ee ee a A SE eR 
; / 
H - 
i = a ic 
: ] . 
2 
{ i 
' 
| 1% 
ee ° —_- 1 
{ aaa Woo if ce, 
| - Mapid 
‘ t 
, —_ = i -- 
~ > 4 i \ 


| | 
eS ee ae | 
. , 
L\ ' 
ae | 
TOW Gree BA <= +446 ~ = ee 
it 
: j 
_ = 4 
ee aad 
ks 
j i - 
| ae tea 
( id ’ 
) 
; alien a a a ors S641 he —_ 
‘ Tt ; 
a < 
{ 
L : \j a soc 


Fi 
a pentane einer eee tks | deter ee “ - , - 
. 


~ f a 
me a ' i 
6 
; _ 
~\. 
° a Py 
3 @ - As, C 


nohiyloeei@ iekbeg : [49 sangeet’ 
; | oT 
Wy . | ss 


i : ‘ i ' J we 


18 


D, 9 D 3C. oC. 
cher tar) hye De Verse 22215) 
Fac) = To 
C.(r,,t> 0) = Cc 2.2(6) 
C. (r,t = 0) = 0 


The instantaneous radial solvent concentration is determined 


analytically using the method of Laplace transform in Appendix A. 


1E 
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: 2.2(7) 
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fret) ae ) = € - *in(iz) | 
Pp = 


The solvent mass within the sphere is obtained by spatially 


integrating the radial concentration. 


ae 2 c A 
Pe Cr dr = C* [oreo + i #COFCE-Dar| 
2.2(8) 
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Substitution of solvent mass into the conservation equation 
2.2(4) yields the radial dissolution integro-differential equation for 
the sphere radius, ony as a function of the macroscopic liquid phase 


solvent concentration, C. 
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9 Ce 2 ; 
ce Be (TRH) +f £ COP (Eta f = tp [ck + acc-cr) | 272(9) 


2.3. Convective-Diffusion with Adsorption 

All macroscopic properties are assumed uniform in the vertical 
direction such that the rectangular lamina collapses into a degenerate 
two-dimensional problem. 

Solvent with constant density, Ps» is transported in the liquid 
phase by both convection and diffusion. The convective solvent mass flux, 


— 


Wo is due to gross fluid movement, 


Me p Cu page OD) 


with the volumetric fluid flux density, u, as the apparent Darcy 
velocity. 

In a flowing fluid, molecular diffusion in a porous medium is 
greatly enhanced by directional velocity dispersion. Dispersive 
behaviour is due to the interstitial inertia of diffusing species along 
tortuous capillary paths. For a point solvent source in a static fluid 
and homogeneous medium, Fickian diffusion results in the radial 
propagation of a spherical concentration front. Superposition of a 
velocity field deforms this spherical shell into an axisymmetric surface 
skewed towards the direction of flow, as illustrated in Figure 2.3-1. 

The principal axis of dispersion coincides with the flow direction 
at every point in the flow field. An instantaneous local velocity frame 
F {xi} can be defined by a longitudinal basis vector, ey» which 


coincides with the flow direction and a transverse basis vector, e}, 
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Figure 2.3-1 : Directional Dispersion 
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which is normal to the flow. 


|u| 273(2) 


In this frame, Fick's first law relates the apparent diffusion 
solvent mass flux, Ww , to the concentration gradient by a diagonalized 
tensor of rank two, pt » whose consecutive elements correspond to the 
longitudinal and transverse diffusion coefficients, respectively. 
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2.05) 
For moderate flow rates, Perkins and Johnston?” find that the 

apparent longitudinal and transverse diffusion coefficients vary with 

the first power of the velocity magnitude. These apparent coefficients 


are determined in terms of total cross sectional areas for cgs units. 


D 
1 
D =p tz |u| 
e 250A) 
Do Fe 
D. = — + .01570d_|u| 
£ F. P 


The constant liquid phase molecular diffusion coefficient, Do» 
is adjusted by the formation resistivity factor, F_, to account LOT 
the tortuosity of the capillary network. The average particle diameter, 
qd? is the grain size and o is a measure of packing inhomogeneity. 

The diffusion mass flux is taveriant © of inertial reference 


frame and can equivalently be expressed in terms of the Cartesian 
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coordinate system. 
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= BESS) 


The similarity transformation involves a rotation”” through a 


u 
oS rant (22] 2 23(6) 
ia 


which determines the matrix of the cartesian tensor. 


velocity angle, 6, 


2 te , 
D D Dd, cos 6+D, Siig (D, -D.Jcos 6 sin 6 


= fi D 2.507) 
D D (D, -D,Jcos 6 sin 0 D, sin o+ De cos 0 


Let ts be a directed surface element on the closed surface,o, 
=~, 
do = ndo PARES Ree. 


where n is an outward drawn unit vector normal to the surface element, 
do. The differential rate of solvent flow through this surface 
element is due tothe combined effects of convective and diffusive trans- 
port mechanisms. The net rate of solvent outflow through the closed 


surface,o, is the integral sum over all elements. 
Ww tw.)+do =o. € ele pli fGi pe WIR MF 2.3(9) 
me Nae oo =e ict ax, | “i 
Oo 


Conservation of solvent mass on the macroscopic scale equates 


the rate of solvent depletion from the combined solute-solution system 
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over the volume element to the net rate of solvent outflow through the 


bounding surface. 


3 ee aC a hes 
= f $ SCdv + ff kA a) Ds; ee Us je. ida 2.3(10) 
AV R. fo) J 


The total solvent content in the solute has contributions from 


all spheres in the volume element, 


ne 
® — 
ff cv. = ff (Nn. f P c_r‘dr) dv = & 6 (1-S)C_dv 2.3(11) 
Rogen é (\, Bee Av ™ % 


and can equivalently be expressed in terms of a macroscopic volume- 


averaged solid-phase solvent concentration, C. 


Fs 3 *D 2 
Coes eal Cor dr 25 E12) 
_ 0 


The volume element is stationary, allowing an interchange in the 


order of integration and differentiation. 
Hy po" (aes (1-Sie_ fav = 1 D.. = -"u.cle,-ndo 2.3(13) 
he ee ot f Ss és ij ue i i 


The rate of solvent accumulation in the solute is related to the 
interphase solvent transfer by the dissolution equation 222(4) 5 such 
that the integral form of the convective-diffusion with adsorption 


equation is expressed in terms of the dependent variables. 
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A 9C An ees 
= ate aes oC] e, * ndo 2.5(14) 
0 J 
Application of the Divergence theorem, 


t 
Dae, 


R) aC 
- ox, Ps ox; - uc |p aw = 0 2.5045) 


4 3\ 9C 
yn { (‘, “ a N ae + 4a x0 [a(C-C*) ] 


determines the differential form of the convective-diffusion with 


adsorption equation for all macroscopic points in R. 
3 dC y 40 3 \ 9C Z Ma 
Ox, Ee Boe | = (+.-4 ‘) pet Anh el Ue Cx.)s| 2 ivi.) 


2.4. Darcy-Continuity 

Convection currents in the streaming solution are established 
by local pressure gradients within the porous region. The apparent 
velocity vector, U4 is the rate of gross fluid movement through a unit 
area normal to the direction of flow. With the assumption of sufficient 
capillarity and flow rates within the laminar regime, the velocity is 


related to the potential gradient by Darcy's law. 


See aes Mee 2.4(1) 


Density variations are assumed small, such that the force 
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1 4 : ; : 
potential Q per unit volume, ¢, is the combined hydrodynamic pressure, 


p, and the hydrostatic gravity head, pgh. 
® = p + pgh 20402) 


The fluid density, p, is determined by the concentration weighted 


average of the component densities, 
aa (pee Jeet 2.4 (3) 


and g is the acceleration due to gravity. 

The permeability tensor of rank two, ie 1s the conductance of 
the porous medium to gross fluid movement and is intrinsically dependent 
on the geometry of the local capillary networks. These networks give 
preferential directions to gross fluid movement. For convenience, the 
principal axes of the permeability trend are chosen to coincide with the 
coordinate axes. 

Dissolution of the solid phase will create a porosity increase 
resulting in capillary enlargement and enhanced flow conductance. 
Similarity of capillary enlargement is assumed for the packed bed and 
the component permeabilities related analytically to the porosity by 
the Kozeny-Carman equation. * 


D4 Mes 
i it a 2.4(4) 


F489 98 c1-0)° 


The permeabilities, k, , are measured at the reference porosity, bo - 
0 
The fluid viscosity, np, is the fluid resistance to flow due to 
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internal cohesive forces. The viscosity is assumed to be independent 
of fluid inertia for the flow rates of the laminar regime and dependent 
only on molecular interaction. 
For the solvent-bitumen system, the viscosity is related to the 
liquidity determined by the solvent concentration, component viscosities 
42 


and densities using Cragoe's Method. The method defines a liquidity, 


L, based on the viscosity measured in poise. 


1000 2 
LQ) = aaa 2.4(5) 
gn yp - gn(5x10 =) 
The fluid liquidity is calculated by the mass-fractional weighted 
average of component liquidities, 
= Bs 2.4(6 
L(u) = MLQ.) + (1-M)LG)) (6) 
with M. as the solvent mass fraction. 
Ps 
Mile -—=-C 24 C7.) 
Su 8/80 


The fluid viscosity is functionally dependent on the concentration, 


with parametric dependence on component properties. 
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Although the individual fluid components are incompressible, there 
can be mass accumulation in the volume element as a result of porosity 
changes from dissolution, and density variations due to mixing. 

DRecoluti on of solid matter leads to the creation of fluid as a result 

of nonreversible phase changes. Conservation of fluid mass dictates the 
rate of mass accumulation over the volume element is a direct consequence 
of the mass inflow through the bounding surface, combined with the mass 
created from the internal fluid-generating source. The dissolution 
Mass-source creates fluid with ambient density approximated by the local 


macroscopic average, at the rate of solute dissolution. 


<— ooddv = - ff pu - ndo ts ff fe) oP 2.4(10) 


The source effectively cancels fluid accumulation due to porosity 
changes, leading to the integral form of the continuity equation. 
ff ¢ 22 av = - ffou- ndo 2.4(11) 
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Applying the Divergence Theorem, 


3 dp 
Ff Mayet Or) dys O 2.4(12) 
AV ee 1 ot | 
determines the differential form of the continuity equation for all 


macroscopic points in R. 


Ay, beU) = >> a 2.4(13) 


The combined Darcy-continuity equation, 


ok. 
Ay i leuteaste, L, At 3 \ 90 
an Go (erent = ne Nr )z 2.4(14) 


OX. u OX, 
together with the convective-diffusion with adsorption equation 
2.5(16) and the radial dissolution equation 2.2(9) ‘constitute the 
coupled system of nonlinear equations for the state vector (p5C5r.,) 
whose closed form solution determines the evolution of the leaching 


process. 


2.5. Boundary Conditions 

A thin surface layer, ) , of thickness 6 separates the lamina 
R from all external points. External preperties are generally different 
from those just within the porous interior, and solvent transmission 
is regulated across the surface layer. The interior of R is void of 
external fluid sources, such that the only introduction of fluid mass 


into the porous lamina is through the bounding perimeter. 
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There is no accumulation in the surface layer and the solvent 
influx from an external reservoir (+) balances the influx into the 
porous interior (-). Mass transport across the surface layer is 


normal to the bounding perimeter and all tangential fluxes vanish over 


) 


oC a oC ; a a 
fie- P45 | = fey, cm ; e. n=0, P(x) « ) 2554 ))) 
dj o1eCry I1C=) 

As a consequence of the normal velocity constraint, the off- 


diagonal elements of the cartesian diffusion tensor coefficient vanish, 


restricting cross diffusion to the interior. 


= 109 ea) eRe 5 2.5(2) 


Three boundary types are considered for the two-dimensional 
rectangular lamina of dimension (L,W,Az). The set of point and/or 
line sources, {} ont constitute all fluid inlets on the perimeter. 


The set of sinks, » constitute all outlets. The complementary 


ene 
set of points, neither inlet nor outlet, are impermeable is FLOW! * 
The normal velocity components vanish at all impermeable points. 

Both concentration and pressure are continuous. Solvent transport is 
restricted. No solvent is allowed to diffuse in, or out, through any 
impermeable point. The normal components of the potential and 
concentration gradients vanish, with Neumann boundary conditions to 
satisfy the imposed impermeable constraint. 

2. Be -0; 8-8 = 0.7) € fy prow 2.5(3) 
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Fluid is removed from the porous region at an outlet pressure, 
Pour: The outlet channel is designed to remove solvent with the outlet 
stream, and filter any solvent from diffusing back into the interior. 
The concentration is continuous over the outlet channel and the fluid 
velocity uniform. The normal component of the concentration gradient 
vanishes, such that solvent passes unidirectionally from the rectangular 
region to an isolated outlet reservoir sink. The fluid front is 


annihilated at the isolated outlet singularity, resulting in an external 


discontinuity within the velocity field. 


SS ee ae 
i 


(+) 2.5(4} 
ae n= 0 PR) € one 


Solvent of concentration CoN is injected from an external 
source reservoir into the porous region through all inlet points at an 


inlet pressure, P Concentration gradients are established at the 


IN° 
inlet surface as a result of the discontinuities in diffusive properties 
between the external and internal boundaries. Solvent diffusion from 
the external reservoir is regulated by selectively passing solvent 
unidirectionally into the interior. In this manner, the source 
reservoir remains isolated and is unaffected by the processes occurring 
within the porous region. Regulation is achieved by allowing solvent 
mass to accumulate on the external inlet surface at a rate proportional 
to the solvent diffusing into the interior. 
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The skin factor, pier is the constant of proportionality which 
determines the diffusion conductivity of the inlet surface. A fluid 
front is created at the isolated inlet singularity causing an external 
discontinuity in both velocity and concentration. The fluid velocity 
is assumed uniform within the inlet channel, and conservation of mass 


equates external to internal fluxes, 


2 fu,- Ds. ma 2.5(6) 
Lace) 
to give the inlet boundary conditions. 
-u. 
1 (=) (1-a.)D, ry : : ee 
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As the skin factor approaches unity, the inlet gives no 
resistance to solvent diffusion and the internal concentration 
approaches that of the external source. When the skin factor vanishes, 
the inlet becomes impermeable to diffusive transport and fluid of 


concentration CN is introduced into the interior by convection only. 
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CHAPTER THREE 


SIMULATOR DEVELOPMENT 


An overall numerical approach to the nonlinear system of integral 
and differential equations for solvent leaching of bitumen spheres is 
presented. The numerical solution for the state vector of fluid 
pressure, fluid phase solvent concentration and the fluid phase 
saturation using the method of successive substitution is outlined. 

Central-differences on a regular mesh centered grid with a 
Crank Nicolson ae discretization determines the finite difference 
analogs for the integral conservation laws. The time step is based on 
the incremental fractional pore volume of fluid injected at a constant 
flow rate, over each discrete injection interval. For line sources, 
the Newton-Raphson method iterates the inlet pressure necessary to 
Maintain the designated flow rates. A mass balance check is used to 
monitor the differential and cumulative errors that may arise in the 


approximate solution to the overall system. 


5.1...Central Differences 

The rectangular region R of dimension (L,W,4z) is partitioned 
into a grid of N, and th equidistant intervals in the x- and 
y- directions, respectively (cf. Figure 3.1-1). Each discrete mesh 


point Pai has the coordinates (x; 5¥5)> 
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with the respective equidistant intervals, Ax and Ay. 


Ax = 
Ke a) 


Ay = 


Each mesh point Ps. has an associated block of volume AV, E 
J ij 
whose bounding surface is determined by the vertical midplanes of height 


Az lying between adjacent mesh points. The midpoints are fixed by 


their spatial arithmetic means. 
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Let i* and j* denote the corresponding faces of those 
volume elements which map onto the boundary perimeter, ) 
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with cross sectional surface areas AY ait and volume AV), 


Leet 
. ; ean 
an is 1 the 


iS i r 
(Syne 
’ j 
ae 
eR. ic 4 nf 
} . : i Py 
i I M ii | | es 
ih anton: ie bnabidoten ab eed ae “aston ian 
i rg <2 be 


treisd %s eongh ipbie baer rhe seal saitneaiall as ‘sgctvue garth 
¥d boxtt ots saat onto ont Jehuton dena! angonibe re 


| if 7 


oy : | eck aati " | 
| | | i 
= s f I f i if we yd) } 4 
Sai ee a ie: 
cE ) ' - i ; Bet 4 i 
F ' 5 ovate 
- ae ' %, ‘ 3 Rad , F A is ay 
| ee ah es me, | ‘ ws ; 
a7 Bas ' . ee 5 : “yg 7 


asoda 3d 26 2 antiaeet ods stone t ba! 4 Gore: 
; i yee ri Pome, 7 
d . pennies +; : 


Me ei a 


35 


Ay 7 AY 582 
J 
ay, = AX Az 5. (6) 
at 
AV, = AX; AY ;Az 
1J 


Consider the nine-point difference star centered about mesh 
point Hae [ef. Fipure 3277-2]. YLet f(x,y be any function ‘which ‘can 
be expressed by a Taylor series in the neighbourhood of ag By 
expanding f about the adjacent midpoints, the central-difference 
approximations to the central first-order derivatives are obtained with 


second-order accuracy. 
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Second-order derivatives are approximated to second-order 
accuracy with the successive application of the first-order difference. 
The adjacent mesh points are used in place of the midpoints to evaluate 


inner first-order differences to approximate the mixed partials. 
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Figure 3-1-2 : The Nine-Point Difference-Star. 
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The mixed partials arising from the cross terms of the diffusion 
tensor product are small, but nevertheless included in the analysis. 
Both first- and second-order pure partials are evaluated using only 
those points on the five-point sub-star, and thereby categorizing this 


central difference scheme as a five-point method. 


».2. Numerical Dispersion 

The component convective fluxes can be expanded to second order 
accuracy at the bounding faces of the volume element. The spatial 
variation in solvent flow over the volume element approximates the 


convective derivative. 
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All dependent variables must be expressed in terms of the mesh 
points, and to the same order of accuracy, the mid-concentrations are 


replaced by their spatial averages. 
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A clever regrouping of the central difference approximation 


isolates an upstream weighted component. 
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The nea dispersion counterpart is termed accordingly for 
its similarity in form to the central difference approximation of the 
diffusion derivative. Chaudhari!” obtains this numerical dispersion 
term directly from a second-order correct upstream-weighted 
discretization. 

Even though central differences and second-order correct upstream- 
weighting give identical discretizations, the latter form will be 
adopted to accentuate the smearing of truncation errors that may arise 
from first-order correct methods, and quantify these errors with a 
physically analogous dispersive mechanism. One should keep in mind 
that second-order correct methods also have truncation error, and the 
weighted combinations of these higher order derivatives will produce 


numerical distortions to some degree. 


3.5. Time Integration 


The time domain is partitioned into sufficiently many irregular 


intervals to determine the required leaching history of the miscible 
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flood. Each time step, At, with respective upper and lower bounds 


eal and oa is based onthe incremental fractional pore volume, 
APV, of fluid injected at a constant injection rate, GIN: 
oy eWdz 
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Consider a nonlinear differential equation with the form, 


Spud 
pee GES) 3£3(2) 


where y,f and re are any functions expressible by a Taylor series 
about a discrete point in time. Numerical integration may be performed 
using a number of techniques, ranging from high-order Runga Kutta 
methods demanding many points for Lagrangian interpolation, to Riemann 
sums using two-point rules. Both midpoint and trapezoid rules give the 
same order of accuracy as central differences, and require a minimal 
amount of computation to evaluate. 

Integration over the time step using the trapezoid rule for the 


Tight hand side of 3.3{2) ‘gives: 
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In comparison, the midpoint rule preserves the nonlinearity, 
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but requires interpolation to evaluate mid-time values. To the same 


order of accuracy, linear interpolation of the dependent variable relates 
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the mid-time value to its neighbouring nodes. 
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This midpoint approach assumes that only the dependent variables 
are piecewise linear, whereas the trapezoid rule averages the non- 
linearity. The author feels that the midpoint rule tracks the non- 
linearity more effectively while requiring less work in evaluating a 
single mid-time coefficient. 

The simulator in discussion adopts this midpoint weighting for 
the temporal discretization of the governing nonlinear system of 
leaching equations. 

The integro-differential radial dissolution equation [2.2(9)] 


LSudiscretized over the time steprto give: 
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The convective-diffusion with adsorption equation [2.3(16)] in 


compact tensor notation is temporally discretized to yield: 
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Midtime pressures are implicitly evaluated in order to directly 
calculate the midtime convective-diffusion coefficients. The 
temporally discretized Darcy-continuity equation 2.4(14) has the 


form: 
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53.4. The Associated Difference Equation for Particle Radius 
Time levels have been established for discretization of the 
radial dissolution equation in 3.3(6). The convolution integral is 
approximated by a Reimann sum using the extended mid-point rule. 
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‘The particle radius is the volume average over the volume element 


AV, and may vary spatially from element to element. The discrete 


ij 
version of the radial dissolution equation for each element is: 
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Isolating for the new particle radius, a » gives the associated 


difference equation for radial dissolution, with the standard form: 
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Iterative techniques must be used to solve the highly nonlinear 
algebraic equation for the particle radius. The summations of F 
require many terms for convergence. Caution must be excersized as the 
dependent variables approach zero, for the equation becomes singular. 
Dissolution is nonreversible and must terminate at the zero radius, 
and hence the integro-differential equation is valid only while solid 
material exists. Once dissolution is complete, a zero radius must 


prevail. 


T rtO 3.4(6) 


3.5. The Associated Difference Equation for Pressure 

The specified boundary conditions will fix an inlet pressure on 
the perimeter to produce the desired injection rate. Flow rates are 
monitored at the block boundaries over each element by the integral 
form of the Darcy-continuity equation. The mass flow transmissibility 
between adjacent mesh points determines inter-block fluid transmission 
across common elemental faces. With the time levels established in 
equation 3.3(8), the central-differenced Darcy-continuity equation has 


the integral form: 
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where both fluid mobility and grid dimension determine the flow 


transmissibility. 
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The associated difference equation for pressure at each mesh 


point has the standard form: 
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Special consideration is given to the mesh points which map onto 
the bounding perimeter, ) - If no inlet fluid sources, or outlet fluid 
sinks, are specified at the perimeter block, then the face which maps 
onto the perimeter is impermeable to flow. Reflection conditions, using 
fictitious points outside the physical region, are constructed from the 
finite difference analogs of the differential impermeable boundary 


conditions. 
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Identical algebraic expressions are obtained for the integral 
difference form by constraining the corresponding fluid mass trans- 
missibilities, and consequently the associated difference coefficients, 


to vanish. 
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If a fluid outlet is specified, the corresponding point and 


subsequently the outlet block is maintained at the outlet pressure. 
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except the central ones to vanish. 
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The discrete Newton-Raphson method is used to iterate midtime 
inlet pressures required to sustain constant injection rates for line 
sources. A line source is any succession of consecutive inlet mesh 
points, with a minimum restricted to two. The inlet faces of these 
consecutive blocks fuse to allow fluid injection over their combined 
Hengih. inlet pressure as constant... adhe, fluid flux. need not be 
uniform over the line source, however the piecewise sum of individual 
rates is required to total the specified injection rate. Pressure 
iterates are obtained from an initial guess and flow rates calculated. 
The change in pressure with respect to the change in flow rate is 
approximated by a backward difference, and the required inlet pressure 


converges within a few iterations. 
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with k as the Newton-Raphson iteration index. 
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No such iteration is necessary for an isolated point source. The 
perimeter faces of the point inlet become impermeable and a mass source 
term equal to the fluid injected at the specified rate is included on 
the right hand side of the Darcy-continuity equation, 3.5(1). 

Discretization has quantized the dependent variables both spatially 
and temporally, such that the step discontinuities at the block boundaries 
and mid-time levels have been linearly averaged out. Nonlinearities 
stem from the implicit dependence of the mobility on both fluid 
concentration and saturation. Mid-time mobilities evaluated at the 
midpoints are subject to both spatial and temporal discontinuity. 
Averaging techniques must be used to weight implicit dependence to 
existing mesh points and time nodes. Such techniques as upstream or 
downstream weighting, and arithmetic or harmonic means, evaluated ahead 
or back in time, may be considered. Each will produce different results 
for different degrees of nonlinearity, and one must numerically 


experiment to obtain the most suitable weighting. 


3.6. Associated Difference Equation for Concentration 

The pressure field, implicitly evaluated at the mid-time, has a 
deliberate phase lag of half a time step in order to directly calculate 
the mid-time velocities. The scalar potential field is well defined at 
these mid-time nodes and the normal velocity components at the bounding 
faces of the volume element are readily obtained using the Cates 


differenced Darcy's law. 
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Even though the normal velocity components are defined at the 
bounding faces, the tangential components may be discontinuous as a 
result of spatial quantization. To avoid slippage along these 
numerical boundaries, the tangential components must be single-valued. 
Both velocity components are defined at the mesh points, allowing a 
Simple arithmetic mean to relate the central components to those at the 


peripheral faces. 
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Cartesian diffusion coefficients are implicitly dependent on 
both velocity angle and magnitude, introducing further nonlinearity 
into the system. As with the mobilities, the diffusivities at the 
bounding faces must be spatially averaged to weight implicit dependence 


to existing mesh points. 
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The midpoint velocity angles and magnitudes can be evaluated 
using similar upstream, downstream or arithmetic spatial weightings, 
as mentioned previously. 

.-Analogous to the mass flow transmissibility, the solvent 
diffusivity controls the solvent transmission by normal Fickian 


diffusion. 
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The numerical dispersivities follow from the central-differenced 


convection term, 
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with the convection transmissibilities as the interblock fluid flow 


rates. 
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Cross dispersion, as a result of the off-diagonal terms of the 
diffusion tensor product, determines the transverse solvent diffusion 


tangential to the bounding elemental faces. 
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With the time levels established in equation 3.3(7), the 
central-differenced convective-diffusion with adsorption equation has 


the integral form: 
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The central accumulation term is spatially weighted over the five- 
point difference-star in order to cancel a portion of the convective 


truncation error. 
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The truncation cancellation procedure, effective when convection 


dominates, determines the optimal weights (w, Wo) which minimize the 
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combined convective and weighted accumulation truncation errors. The 
combined error, consisting of third-order mixed partials, is 
transformed into the velocity frame, where the derived weights suppress 
a dominant longitudinal error term. Appendix B gives a detailed 
derivation in which these weights are analytically determined for a 


direct solution scheme. 
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for an alternating direction implicit (ADI) procedure. 
The truncation cancellation procedure is limited to a maximum 


time step determined from a stability analysis in Appendix B. 
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convective dominat problems associated with miscible displacement. By 
setting both weights to zero, the proposed discretization collapses to 
a second-order correct central-difference scheme, applicable to the more 
general class of leaching problems. 

The resulting associated difference equation for concentration 
has the standard form: 
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The cross dispersion terms are included in the right-hand side 
of the associated difference equation and their values are subsequently 
iterated with all other time-dependent coefficients. 

As established, the perimeter is impermeable to solvent unless 
an inlet or outlet is specified. The normal components of both 
velocity and the concentration gradient must vanish at all impermeable 
faces of the perimeter blocks. Reflection conditions yield the 
superposition of a negative flow field, in which a reverse-running image 


concentration-wave accounts for the total reflection of solvent at the 


no-flow boundaries of the finite porous region [cf. Figure 3.6-1]. 
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Figure 3.6-1 : Impermeable Boundary 
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Identical algebraic expressions are obtained by constraining 


all transmissibilities to vanish at these no-transport boundaries. 
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Tangential fluxes exist on the interior of the impermeable faces 
and only the normal velocity components vanish at these points. 
Solvent enters and leaves the interior normal to the bounding surface, 
thereby the tangential fluxes vanish at all inlet and outlet perimeter 
mesh points. As a consequence of the velocity angles at the perimeter, 
the off-diagonal cartesian diffusion coefficients, and subsequently the 


cross dispersion terms, vanish at all perimeter faces. 
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Fluid is introduced into, and removed from the porous region 
at a uniform rate and hence both inlet and outlet velocity gradients 
vanish. Unlike the no-flow boundaries encountered at the impermeable 
points, the inlet/outlet points involve the partial transmission of 
solvent. Superposition of a negative velocity field results in a 
Singular discontinuity, in which fluid fronts are either created (inlet 
source) or annihilated (outlet sink) to preserve the material balance. 
These source-sink terms must account for the fluid distributed over both 
the physical and image counterparts. 

Solvent is removed at all outlet points with the fluid stream by 
convection only and thereby the normal component of the concentration 
gradients vanish. The resulting finite concentration wave is obtained 
by the superposition of a diffusion reflection defined in the physical 
domain of the porous region [cf. Figure 3.6-2]. Reflection conditions, 
including a distributed solvent sink density with a superimposed 


negative flow field, resolves the fluid singularity at the outlet point. 
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Sink Boundary 
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The corresponding outlet transmission conditions allow solvent 
to accumulate at the outlet face, as if impermeable to flow [3.6(17)]. 
The inclusion of a solvent sink equivalent to the rate of solvent 
outflow provides the necessary fluid removal where diffusion out of 
the porous region is fully damped. The sink, distributed over the 


outlet volume element, 
n+43 _ nt+%3 n+4y\ nts, p 
See ina (omen Uy )er or tue our? 3.6(20) 


is included on the right-hand side of the associated difference 
equation. The factor "two'' in the sink density of 3.6(19), arises 
from the half- and quarter-blocks on the perimeter. 

Fluid, with solvent concentration Cine is injected into the 
interior through all inlet points. The skin effect, includes an 
accumulation term distributed over the inlet surface, proportional 
to the solvent diffusing into the interior. The effective solvent 
source depends on the diffusivity of the inlet surface, which is quanti- 


fied by the skin factor. Reflection conditions, with a superimposed 
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negative flow field and an effective solvent source density, resolves 

the fluid singularity at these inlet points [cf. Figure 3.6-3]. A 
concentration discontinuity may exist at the inlet surface, where the 
external concentration is fixed at CoN? and a non-zero gradient just 
within the interior. Reflection conditions which infer the impermeability 
of the external surface, in conjunction with the skin accumulation tern, 
satisfy the non-zero gradient internal boundary condition. In this 
manner, the external source is an isolated singularity and gradients may 
exist along an injection path, be it the dimensions of a well bore or 


Ounthe microscopic scale of the grain size. 
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Figure 3.6-3 : Inlet Source Boundary 
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The transmission conditions for the inlet are those of the 
impermeable boundary [3.6(17)] with an effective solvent source 
allowing for the partial diffusive transmission from an external 
reservoir. The effective source, distributed over the inlet volume 


element, 
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is included on the right-hand side of the associated difference 
equation. 

The concept of "skin effect" provides mathematical flexibility, 
in which the Dirichlet condition is a limiting case. Diffusion of 
source solvent in a point configuration, representing an injection well, 
may be negligible for a large field yielding a virtually non-existent 


skin effect. In contrast, a line source representing fluid injection 
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into a linear core, has significant surface area open to diffusion and 
consequently a diffusion source must be included. 

An initial surface concentration discontinuity may exist at all 
inlets if the internal solvent distribution differs from the external 
solvent concentration to be injected. The discontinuity is resolved 


+ 
by designating the initial surface concentration, 0° » as the arithmetic 


F pepo: : 0 
mean Of the initial block ‘concentration, ‘C°”, ‘and ‘that of the 


concentration to be injected. 


oy =5 a as 3.6(23) 


The limiting Dirichlet case requires, 


ees. yee 3.6(24) 
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and the mid-time injected concentration follows deductively. 
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apply these limiting assumptions to a physical problem. The injected 
concentration or is assumed to be constant for all times, however 
; n+s : 
the mid-time injected concentration Cry” remains dependent on the 
F n 5 
current surface concentration, Ga . In this manner, temporal 
stability is maintained as the surface concentration approaches the 


injected concentration. Stiffness of the inlet condition is ensured for 


skin factors approaching unity, however smaller skin factors may generate 
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either asymptotic or oscillatory responses depending on the inlet 


Peclet number. 


S47 s \ Direct’ Solver 
Each state vector has an associated difference equation for each 
point in the mesh centered grid. Both pressure and concentration include 


spatia) difference equations with the standard form, 
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The associated difference system, ADS, is the system of N 
difference equation with N_ scalar unknowns, each corresponding to 


the central mesh value of the dependent variable. 


N= N_N Or ey) 


The natural ordering index, <j, maps a two-dimensional mesh 
point onto a one-dimensional array vectored sequentially along consecutive 


sweeps of the grid [cf. Figure 3.7-1]. 
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The ADS is a vector equation, 


= 3.7(4 
Anat *nat Qrat 1) 


where the non zero elements along each row of the Nx N_ coefficient 


matrix corresponds to the associated difference coefficients about each 


natural mesh point. 
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Anat AL) = ay 

Anat ts t4t1) = Ba a ee ee N 3 uv {mod(éj ,N,) #0} 
Anat 4d ti-1) 7 Dy, 3 TE {af 1} \ {mod (éf N,) #1} 
Anat ts 44ND = Hyp 3 if ra #N,} tie BONN 


Anat Us» NU = Bee Pelee ls tape 


3.705) 


The rectilinear constraints confine the numerical system within 
the spatial bounds of the physical grid. The modulus operator converts 
the global position to a local sweep value where the equivalent natural 
conditions ensure a square matrix. 

Conventional Gaussian elimination techniques require the storage 
of all elements and the execution of numerous row operations to reduce 
the coefficient matrix into row echelon form, before a backward 
substitution can generate the solution. Better resolution requires 
excess mesh points forcing the augmented matrix to grow exponentially. 
Gaussian elimination becomes time consuming and clumsy, with many 
Operations involving zero elements. 

With the structure of the matrix map well established, an 
efficient solution algorithm is constructed. The scheme utilizes 
symbolic mappings to minimize the work necessary to generate a solution, 
while simultaneously reducing the storage by accomodating nonzero 
elements only. 

By reordering the grid along alternate avaconare** icf. Figure 


3.7-2], the coefficient matrix partitions into well defined quadrants. 
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This partitioning nature of D, ordering, allows an L-R decomposition 
to be applied with relative ease. A coupled recursive procedure 
requires minimal work space to unravel the upper and lower triangular 
matrices. By first using a forward substitution to determine an 
intermediate vector, a backward substitution follows to obtain the 
solution. 

The D, map assigns a pointer array which maps each natural 
index to its corresponding D, index. The inverse D, mapping uses 
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a separate pointer array. 
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has a transformed coefficients matrix whose nonzero elements scatter 
into well defined quadrants. 
Each row of the transformed coefficient matrix will have new 


associated coefficients which must abide by the equivalent Dy conditions. 


‘8 


mort heotpmasb Bl “ite ox sebsinics ao to ree 
arckitya semi ‘pelieniton beirquos Ls can ie it 
0B ‘gyn, Tet sh teaqgu os fovarnmay 03, 
6 an teteite os PORTE Fedice ar a a 

att’ ihagiel iy onion mort ins enh: ty 
(erie aise Saale vinta ‘Gitte ASNLAg ® ' 


even gadeggom 2, Mewavat seit Pak gay i hier 


= 


ASP 
i | 
eS ae 
(ss ee (>! pep: . 
" ' 
conta 


co] 
pa 


TI Tene rruietne £2 ea sirrgut ag xiTIEM on loti. bowroRenat? ® ‘aad a 


| a -2inerbewp beniten {ew osm _ 
watt, e¥erd- i biw xine sani aise tanoenee ant 0 wor soe 


diye Pa thaws: 4 29 Jaslaviupahiedy abate Seats as bale ae, besaisoess , 


ae ; : : 1 nt 
- - U ; : ae 7 * | 


68 


Ay Ui ii) = BD,'<Kf>) 


ANNES .D,<D, <Lf> +12) F(D, <is>) ; 
4 


é ait be 
if { mod (D i<éj>.n) #0} 


Mis lee ale 
An, Das D, <4Aj>-1>) D(D, Aj); 


: a a ee 
if {mod (Dy <A{ 7 N,) #'1 \ 


ee Se 
Ty, eb enieae plain 3) H(D, <Af>); 


gine 0, kéj> SN \ 


ul 


Ber a Oe = leer 
uy ESE aay yen) B(D, <A4f>); 
: Sy Os 
Stjo 
er {0, Af Nn} 
3.7(9) 
A specific example for the (5% 3)) -orid-is shown in Figure 3.7-3 
CO 11 lustrate® the partatsomne nature wrethe D, transformation. The 
general form of the transformed coefficient matrix has four submatrices 


[cf. Figure 3.7-4], with the upper left and lower right quadrants 


diagonal. 


An L-R_ decomposition of Ay will yield upper and lower 
4 


triangular matrices, 


A. -= "LR 537416) 


with the general forms illustrated in Figure 3.7-5. The lower 
triangular matrix, L, has diagonal elements of unity. The L-R 


decomposition can be equivalently expressed in terms of the component 


submatrices. 
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Since Lay is the-identity matrix, 


ars Slay 


it follows directly that the upper submatrices of R, are identical to 


the upper submatrices of An J 
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With Ri identical to the diagonal submatrix, Ay. Loy follows 


deductively. 


Pes 2 3.7(14) 


= 
Loy 4 AoA) 


Inverse diagonal matrices pose no difficulties to evaluate and the 


row-wise elements of Loy are determined with relative ease. 
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difference coefficients, and the manner in which they scatter into the 
workspace. The development is based on inspectional pattern recognition, 
with a sequence of tests assigning groups to specific locations. An 
unravelling procedure recursively restructures the workspace, to be 
used in conjunction with both forward and backward substitutions. 

The lower triangular submatrix has diagonal Lo whose elements 


are unity. Intermediate quantities, L* and R*, are defined in terms 


of known submatrices. 
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Furthermore, a maximum of four other combinations are possible 
for the remaining nonzero entries of Ro. These combinations, and a 
Similar set for Lo ,» are listed in Table 3.7-1. Each combination has 
an associated test which generates a value of & for the current index 
Aj. If this value for & is between designated limits, the 
corresponding combination successfully occupies the unique position 
(2,44). No two combinations can occupy the same position. Table 3.7-2 
lists these tests. 
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workspace of Ri and Loe With RO and Le fully evaluated, a 


coupled recursive procedure restructures the workspace to unravel 
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lower triangular matrices, the solution of the ADS, 3.7(7), can be 
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The outlined solution scheme was originally developed for the 
general three-dimensional system of coupled difference systems for 
multi-variable state vectors. The general approach involves associated 
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for this scalar two-dimensional case. 


3.8. Successive Substitution 

A leap-frog type solution procedure, similar to that suggested 
by Peaceman and Rachford, !* is used to update the state variables for 
each successive time node. 

The initial iterates of concentration and particle radius for the 
next successive (new) time node is set to their current (old) values. 
The next iterate for the new particle radius is determined pointwise to 
establish the new fluid saturation. The current iterate for the mid-time 
pressure field is obtained from the direct solution of the pressure 
associated difference system based on the current iterates for both 
concentration and particle radius. The mid-time velocities and 
diffusion coefficients are subsequently used in the concentration 
associated difference system to generate the next iterate for the new 
concentration. A convergence test completes the current iteration cycle 
with a comparison between new and old concentration iterates. The next 
iteration cycle restarts with updated iterates. The last iteration 
cycle of successive substitution is marked by the convergence of the 
concentration profile within the limits of a specified tolerance. 

The nested iteration loop with iteration index m is flow 
charted in Figure 3.8-1 to illustrate the successive leap-frog 
substitution. For line sources, a discrete Newton-Raphson iterates the 
mid-time inlet pressure required to sustain the desired injection rate. 
The Newton-Raphson iteration loop (k) is flow charted in Figure 3.8-2 


using the injection rates calculated from iteration loop (m) . 
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Figure 3.8-2 : Newton-Raphson Iteration Loop (k) 
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At each time node, both differential and cumulative solvent 
balance indicators monitor gross computational errors. The indicators 
are obtained directly from global mass balances executed over discrete 
time steps. 

The total solvent volume, S, in the fluid stream at any time is 


the sum over all blocks. 


S= {¢ SCAV, } So8(1) 
ih a 


The total solvent volume injected into the porous region, Sin? over a 


time step is the sum over all inlet blocks. 


a 
n+5 _ Ss n+35 
aa = At eh ) (Au, Ay) lest aay Cw |} 5.02) 


The total solvent volume produced from the porous region, Sour? over a 


time step is the sum over all outlet blocks. 


n+35 
oa an BSS 
= one a By yee (3) 


The total solvent volume transferred into the solid phase, Scor> over 


each time step is the sum over all blocks. 


N } 
n+ n+, 2 | n+] on nee 3.8(4 
= Cs Hr Oct Mitte 7ete CP -8(4) 
SsoL 7 hv ) (r, p n+ Ss ij 


A differential mass balance over each time step equates solvent 


accumulation in the fluid phase to injection, production and solvent 
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transfer into the solid phase. 
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Similarily, a cumulative mass balance monitors solvent flow over the 


total time elapsed. 
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Both differential and cumulative balances are rearranged and 


expressed in terms of mass balance indicators. 


n+l n+45  n+45 
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For perfect mass balances, both the differential indicator, MB > and the 
cumulative indicator, MB A» must equate to unity. These indicators 
monitor gross computational errors only, such that unit indicators do 
not necessarily ensure the correstness of predicted state vectors. 

The time-advance loop (n) is flow charted in Figure 3.8-3 


to illustrate the overall solution procedure for updating the state 


variables. 
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Figure 3.8-3 : Time Advance Loop (n) 
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CHAPTER FOUR 


MISCIBLE DISPLACEMENT 


The absence of interphase solvent transfer andconsequently a 
constant fluid saturation determines the special case of the miscible 
displacement of an initial bitumen distribution. Under these conditions, 
the solid phase becomes inert and the permeability time independent. 

The analytical solution for the linear problem provides the means 
for testing the accuracy of numerical methods. The degree of convergence 
of the finite difference solution to the exact solution is governed by 
the characteristic ratio of convection to diffusion rates and the ratio 
of a dimensionless time step to the normalized cell length. Stable two- 
dimensional concentration surfaces are predicted using the guidelines 
established from the linear problem. 

Unfavourable displacements generate grid orientation and adverse 
mobility related instabilities. Realistic concentration fronts can be 
predicted using compatible mobility approximations for specific grid 


orientations. 


4.1. Linear Convective Diffusion 
The relevant equations describing miscible phenomena are the 


convective-diffusion and Darcy-continuity equations for the state 
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For a linear homogeneous horizontal system with unit viscosity and 
density ratios, the fluid velocity becomes constant. The miscible 
displacement of an.initially uniform bitumen distribution in this linear 
system is described by the one-dimensional constant-coefficient 
convective-diffusion equation with a Dirichlet condition at the inlet 


(unit skin factor), and a Neumann condition for the outlet. 
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the transformed convective-diffusion equation has the dimensionless form, 
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where the Peclet number, Npe? is the characteristic ratio of convection 
to diffusion rates and the dimensionless time, +, corresponds to the 
fractional pore volume of injected fluid. 

The exact solution to this one-dimensional problem is obtained 
using the method of Laplace transform in Appendix C. The analytical 
solution has the form of an infinite series of error functions where 
the successive terms arise from superimposed diffusion reflections at 
the outlet. The series converges rapidly for the Peclet numbers and 
times of interest, and can be accurately approximated by the first 


three terms. 
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Figure 4.1-1 illustrates the lower limit at which the truncated 
series solution is valid. It is shown that for a unit Peclet number, 
convergence is sluggish and more terms may be necessary at larger 
dimensionless times. 

The successive terms are a consequence of the region being 
finite. The first term is recognized as the solution to the semi- 
infinite analog, with the outlet at infinity. Figure 4.1-2 shows the 
convergence characteristics at higher Peclet numbers. The deviation 
of the finite eolaeion from the semi-infinite case is shown. At high 
Peclet numbers, the deviation is small and for intermediate times, the 
concentration profiles may be represented by the first term only. 

The effect of Peclet number is graphed for the family of con- 
centration profiles at the half-pore volume injection time in Figure 
4.1-3. The curves illustrate the formation of a sharp front at higher 


Peclet numbers. 


4.2. Numerical Response and Tracking 

The second-order correct central-difference scheme (CDS) and 
the truncation cancellation procedure (TCP) have been proposed as 
solution techniques. Both schemes stem from term by term central 
difference discretization. Other schemes, including the operator 
compact implicit!’ method (OCI) and the single-cell high order 
echened: (SCHOS), use constraints based on the highest-order accuracy 
possible on the smallest difference-star, to determine suitable 


weighting coefficients. OCI derives a spatial operator for the combined 


convective and diffusive terms using Taylor series. SCHOS, based on the 
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local power-series solution, is intrinsically multi-dimensional and a 
one-dimensional analog is not available. OCI is outlined in Appendix D 
for the linear problem, such that comparison with an alternative 
approach is possible. 

In a two-dimensional nonlinear problem, variable velocities give 


rise to a range of cell-Peclet numbers, N oe. 


Ne eee 4.2(1) 


The behaviour of numerical methods for the linear problem will yield 
valuable insight into the expected performance for two-dimensional 
problems. The numerical solution to the linear problem using CDS, 
CPT and OCI are compared to the exact series solution for a range of 
Peclet numbers. Time steps are varied. 

Figure 4.2-1 illustrates the damping effect of smaller time steps 
On an inherent source oscillation at low cell-Peclet numbers, Noe = 9-1). 
The oscillation originates from the initial concentration step function 
and is localized about the source at low convection rates. Both CDS 
and TCP track the exact solution adequately, however gross deviations 
are apparent for OCI at larger time steps. It is anticipated this 
deviation will cause OCI methods to predict lower concentrations at 
the stagnant pockets encountered in a two-dimensional displacement. At 
these low Peclet numbers, the CDS method provides the best overall 
accuracy. 

These oscillations disappear at moderate cell-Peclet numbers, 


Noe = 0(1.0), and dimensionless time steps not exceeding the normalized 


3 bas raaofensatS. ikon Wtssseatesns at ee : 

@ xibasqqaA ni bent t te at 130 voldel Lave toa at gotens 1 
*yitanrresie me aetw s9e bxsquo> teat aug ott nt 

ik von 2k a 


ovis esttivtey ania tiete tent laon ‘enotensmit-oed * 
wag , eh bdnail solosdatlea i 


cae — 


rt 


blese 41 iw makdeng waonctl eds Oo? ebortem Leo iremust to tuo 
leqoteuemib-awe tok songmidiveq besoeqxze ods coat adgieni aid 
202 ane meltiong tHenh. ald aT xo lwide Lal seabhe eit oP 
to Suhteh s 40% Aob¥eloe etivee\agexe ods oF tenons? ore | 0 1 ‘See. 
| hate ane eqese enit i 

etesite emis t3ilead to t90R%. gnicnab sas tooauredi ti MoS, orugit — 
Ch. 0> ha ileteianead Retuety tins wat ss wATeEl ine somes Pane 
neat teat gore Hebei eich ipcbies vis wo tomolgize rlsntio 


2083 tad /2990% MOtT>ovEDD wol se sounipe alld suode bestieze! att 


wry 


_ 


enot tetyeb 22098 Teves , wlotaupote aod tehar FIERO ony joors oT ie 
ae 
aida se i cis aa pe a erste: Sarts rogiel #8" tno xo Inoraqga ema | 


1a 


a8 enutiee iinet ‘Tewos arate ot ehodgem 190 saues, Iliw 0, 


4 bs 
tA. nigasoal@etb. fenaiensmth ows & nt bexesmgone eredgeq. samagese sit. 14 
5 


| ae 
listeve teed od aebiverq hodtom “20d ens uli ai20% wot sods “RR 
4 | XX i 109 8 


. - 


erent valoe4-Liso ararsbor “38 Teeqqne ib anoi set ttuea aeodT r -- 


bostiawron ada grifansns ape oe. omit? = bre “t= = at 


e Z gem, 
an} - ; 
ee ¥ 


Concentration 


Concentration 


iz 


.O. 2nd — orde 
re} tee 


1.0 


@ 
oO 
© 
oO 
wt 
© 
is 
oO 
Oo 
0.0 0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9 4.0 
(a) 
iV) 
= 


1.0 


0.6 0.8 


0.4 


0.2 


0.0 


Figure 4.2-1 (a), (b) 
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Figure 4.2-1 : Numerical Solutions - Low Cell-Peclet Numbers 
(a) ,(b) Concentration Profiles 
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cell length. All schemes track the exact solution accurately with TCP 
and OCI errors approaching zero [cf. Figure; 452-215 sities expected 
that all three methods will exhibit a well-behaved, accurate 
representation of the two-dimensional concentration surface for 
displacements simulated at moderate cell-Peclet numbers. 

Gross errors are encountered with the OCI method at cell 
numbers above the invertibility criterion, Noe = 3.464 [cf. Appendix D]. 
Figure 4.2-3 shows the result of overstepping this criterion with the 
onset of gross aie instability. Figure 4.2-4 illustrates the 
existence of an upper bound to the instability interval with the revival 
in OCI marked by the damping of temporal oscillations. A stability 
analysis by Ciment, Leventahl and Weinberg?” predicts such an instability 
interval, above which OCI has an oscillatory convergence to the exact 
Solution. Caution must be exercised in applying OCI methods to two- 
dimensional nonlinear problems, where the critical cell-Peclet interval 
need not be identical with that of the linear problem. 

A convective-dominant oscillation is characteristic of difference 
solutions at higher Peclet numbers. The oscillation originates from 
the source and grows with a propagating front. In Figure 4.2-5, the 
response of the solution methods to high Peclet numbers is illustrated. 
A bounded oscillation persists with the second-order CDS for all 
possible time steps, consequently lower concentrations are predicted at 
the front. Reduction in time step provides temporal damping for OCI 
schemes, giving an adequate representation of the front with only a 
slight ripple. The TCP method was originally derived for convective- 


dominant problems, and in this one-dimensional case gives high accuracy 
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Figure 4.2-2 : Numerical Solutions - Moderate Cell-Peclet Numbers 
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Figure 4.2-3 : 
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Numerical Solutions - High Cell-Peclet Numbers 


103 


eof 


4 
i 
c 
Ot ao 
~*~ © = 
et pe All 
) 
Ste a, oo take : 
°° * a Wis @ oy nee - me “a 
Re cn | ee - 
y ; i " “#99 Wr dnetee 
ey | ee Fai. ) 
a “fg } in m ' 
; ae t : : : ae 
‘ iy 7 A 
ee , « a 
' a 
i ¥ i a , fc] ; ~ > 7 2 : in _ 
; a I ‘ ; of 
; ; Lali na sae eee Seana vo & 
| denna °F" ~ 7: 
7, - eee ge : ae of ‘poses saa 
o£ : F a Kank , F : - 
| jo Th 
i 
= i hl ~, a BN 
:: ¥ iv 
% ; i 
' ne I 
: ; 
| us =f 


esedeut! sas sail xan soins eaten: 


i we 
‘pe tee 


104 


for a unit time-step to cell-dimension ratio. 
These general observations are valid for finer grids. Finer 
grids provide a means to simulate at higher flow Peclet numbers, Noe? 


for a fixed stable cell-Peclet number, N Mesh refinement improves 


Cem 
the accuracy for all three methods, however at the expense of smaller 
time steps. A minimum grid density is required for OCI to function 
properly, whereas no set restriction is imposed for CDS and TCP. 

All methods may exhibit small-time deviation. For convective 
dominant problems, this small-time behaviour results in concentrations 
predicted outside the physically acceptable bounds. The small-time 
effect is a numerical response to the initial discontinuities in the 
discrete representation of the propagating solvent-oil contact surface. 
All methods show similar small-time periods over which the numerical 
solution evolves into the stable profiles which track the exact solution. 
Tampering with small-time evolutionary process by truncating the over- 
shoot will definitely produce erroneous results. 

For this one-dimensional analysis the cell-Peclet number is the 
similarity group which governs convergence behaviour — be it absolute, 
oscillatory or divergent. The ratio of dimensionless time-step to 
normalized cell length governs the relative error from the exact 
solution where the absolute minimum error is set by the grid density. 

A simple rule of thumb for the operating numerical conditions 
necessary to produce confident two-dimensional results can be established. 
The physical problem will produce flow-Peclet numbers for the geometries 
considered. An operating maximum component cell-Peclet number is 


chosen such that the specific method has absolute convergence for the 
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associated linear problem. The grid density required for simulation 
follows from the operating cell-Peclet number which in turn sets the 


upper bound for the dimensionless time step. 


A 5f . LHeuFIVe Spot 

Current simulation problems deal with multiwell flooding patterns 
for optimal recoveries. The five-spot pattern is the basis for many of 
these AcE EY rd eee The quarter section element of symmetry may 
be mathematically represented by a point source and a point sink placed 
on the opposite corners of a square region bounded by an impermeable 
perimeter. 

Spatial variations in the velocity field and the anisotropy 
of dispersion make it virtually impossible for any one scalar group to 
govern displacement similarity. Furthermore, nonlinearities would 
render a general Peclet tensor surface ineffective in scaling the multi- 
dimensional convective-diffusion equation. Since a general similarity 
group is not available, dimensionless subgroups influencing displacement 
must be isolated. 

In order to define a dimensionless ratio of convection to diffusion 
rates, a characteristic length is essential. Normal cell-Peclet numbers 
may be defined pointwise using the grid partition. The relative 
magnitude of individual cell numbers depend on both block size and 
position, while endpoint velocities determine an absolute maximum. 

For a fixed grid partition, the injection rate governs velocity 


maxima. An inlet cell-Peclet number follows readily from a symmetric 


partition. 
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u Ax q 
Nee -( - “205 a 4.3(1) 
IN xx JIN XX 


Flow similarity requires partition independence and an inlet number is 


defined with the physical field dimension as the characteristic scaling 


length. 


q 
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The major factor influencing the shape of the velocity field and 
subsequently the Peclet curvature is the fluid mobility. An incompress- 
ible fluid with constant mobility flows along the streamlines of a 
Laplacian pressure field. For a homogeneous medium, the concentration 
dependent viscosity determines the fluid mobility. The ratio of 
component viscosities in the binary mixture provides a parametric 
measure for the mobility-induced departure from a Laplacian pressure 
field, and its effect on displacement similarity. 

As the grain size approaches zero, dispersion diminishes leaving 
diffusion isotropic and constant. Dispersion effectively increases the 
apparent diffusion proportionately with the local velocity magnitude. 
The Peclet surface responds to dispersion by a global drop in magnitude 
with pronounced reduction about the local maxima. 


A longitudinal Peclet number, 


4.3(3) 


can be decomposed into molecular and dispersive counterparts. 
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The longitudinal number asymptotically saturates with increasing 
velocity magnitude. Higher injection rates elevate global velocity 
magnitudes aes that local contributions from the molecular number 
become negligible, and the constant dispersion counterpart determines 
displacement similarity. 

A laboratory scale model with the physical properties listed in 
Table 4.3-1 has been chosen as a basis for simulation. The dimensionless 
time, t, corresponds to the pore volumes of fluid injected into the 


square region. 


ee ae 4.3(5) 


Figure 4.3-1 illustrates the effects of injection rate on 
displacement similarity in a Laplacian field with isotropic diffusion. 
Each grey-level frame represents the concentration distribution with 
solvent (white) displacing oil (black) at the corresponding injection 
time, t. At low injection rates, diffusion dominates with marked 
radial propagation at both inlet and outlet. Higher injection rates 
bring forth convective dominance with the formation of a sharp front in 
the concentration surface. The radial advance in the early stages is 


succeeded by the delayed breaking of the solvent front to the sink. 
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Table 4.3-1: Miscible Displacement Simulation 
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Figure: 4,3-1 ; Displacement Similarity - Injection Rate Effects 
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Higher convection rates reduce the sweep efficiency with larger oil 
deposits remaining at breakthrough. 

With the injection rate fixed for convective dominance, 
Figures4,5-2eLldustrates the effect of dispersion. Increasing the 
particle diameter instigates the formation of a sharp cusp at the 
leading edge of the front, prior to breakthrough. Further increase 
in injection rate, sharpens the cusp till asymptotic similarity 
prevails [cf. Figure 4.3-3]. 

Solvent injection of .0025 pore volumes over each time step 
was used for the (21x21) partitions. Figure 4.3-4 shows the TCP 
predicted concentration surface for a convective dominant displacement, 
prior to breakthrough. Small-time oscillations were completely damped 


within the first few time steps to give this stable representation. 


434, Grid Orientation and Viscosity Effects 

A diagonal grid represents the quarter section element of symmetry 
with the injector-producer geodesic along the square diagonal, oriented 
at 45° to the grid lines. The physically identical flow system can 
also be represented by a parallel grid, with the grid lines either 
parallel or perpendicular to the injector-producer geodesic. The parallel 
grid corresponds to one half of a five spot pattern, where the grid 
boundaries map onto the symmetry planes defined by the geodesics 
[cf. Figure 4.4-1]. 

The nonuniqueness of finite difference solutions for different 
grid orientations was first demonstrated by Todd et Py fisos Finer grid 


definition does not necessarily alleviate the grid orientation problem, 
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Figure 4.3-2 : Displacement Similarity - Dispersion Effects 
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with diagonal and parallel simulations predicting two different 
displacement histories. Adverse mobility ratios aggravate the 
deviation, which is most apparent in five-point finite difference 
schemes. 

Numerous attempts have been made to eliminate the grid 
orientation effect, none of which is entirely satisfactory.°° All five- 
point schemes experience orientation problems at higher viscosity 
ratios, independent of the mobility approximation used. 

Aside from the nonuniqueness, frontal instabilities may develop 
in either orientation. An early indentation of the leading displacement 
edge may evolve into numerical fingers [cf. Figure 4.4-2]. Mild 
numerical fingering produces an inverted cusp which may disappear if 
diffusion is sufficiently large to equalize transverse gradients. As 
convective forces dominate, numerical fingers develop about the 
symmetry plane, leaving unswept oil along the diagonal geodesic. These 
anomalies are not exclusive to five-point schemes, and similar finger- 
like features have been reported for nine-point methods.>" 

When diffusion is small, different modes of numerical fingering 
can be observed [cf. Figure 4.4-3]. The solvent preferentially moves 
along the confining walls, with multiple secondary fingers at higher 
viscosity ratios, bypassing the central oil in place. 

The unstable phenomenon was generated with arithmetic approxi- 
mations for mean mobilities and diffusivities on a diagonal grid. 
Similar effects can be observed on a parallel grid with breakthroughs 
predicted earlier for a wider stability range. Parallel grids are more 


attractive for moderately adverse viscosity ratios, however instabilities 
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Figure 4.4-2 ; Numerical Fingering 
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Figure 4.4-3 : Modes of Numerical Fingering 
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render arithmetic mobilities ineffective for realistically predicting 
grossly unfavourable displacements. 

Mobility-diffusivity approximations, compatible with the grid 
orientation, are sought for the realistic predication of displacement 
fronts over a range of viscosity ratios, extending into the tar sand 
regime. The "realistic" criterion excludes any displacement whose 
leading edge Beye ices fromthe fastest streamline, defined by the 
Minimum energy geodesic connecting injector to producer. Simply, the 
occurrence of any frontal instability resulting in numerical fingering 
is rejected. 

Mathematically, midpoint mobilities evaluated at the midtimes 
are allowed to take on any value within a bounded uncertainty interval 
(Ax,At) [cf. Figure 474-4]. Similarly, midpoint diffusivities. lie 
within a spatial uncertainty interval fixed at the mid-time. Suitable 
weightings are required to relate mid-values to existing grid points 
in space-time. Many options are available, with specific weightings 
labelled in Figure 4.4-4. 

Numerical experimentation isolated the downstream, back in time 
approximation with diagonal grid orientation as the only compatible 
configuration for realistic displacements, over the desired range of 
viscosity ratio. Any other weighting, or change in grid orientation, 
resulted in some anomaly for the displacement conditions investigated. 
The compatible approximation predicts stable fronts with channel 
formation at extreme ratios [cf. Figure 4.4-5]. Realistic displacements 


show diminishing sweep efficiencies, with advanced breakthroughs, at 


higher viscosity ratios. 
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Figure 4.4-4 : Mobility-Diffusivity Approximations 
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Table 4.4-3: Miscible Displacement Simulation 
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Figure 4.4-5  ; Displacement Similarity - Viscosity Effects 
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Table 4.4-4 : Displacement Similarity - Viscosity Effects 
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Figure 4.4-6 illustrates realistic viscous fingers at lower 
diffusion rates. Fingers develop early in the displacement with the 
central finger outgrowing the others. Once the finger develops, Lt 
does not grow to a point but forms a rather blunt front. The same 
statements are quoted by Habermann” for an experimental study conducted 
in 1960. Although the simulation model-characteristics differ, the 
predicted results have definite similarities to specific experimental 
runs at considerably lower mobility ratios. Of these similarities, the 
most prominent is the early separation of the fingers from the confining 
walls. 

Even though stable fronts are predicted, the existence of a single 
approximation compatible with a specific orientation does not guarantee 
the validity of simulated results. One must keep in Mind, these. are 
approximate solutions only, and simulated results must be scrutinized 


thoroughly, before acceptance. 
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CHAPTER FIVE 


LEACHING 


Single-cell dissolution provides insight into the phase 
behaviour of the closed solid-liquid system. Two mechanisms govern 
dissolution kinetics — solvent adsorption and solvent absorption. 
Adsorption is the surface phenomenon responsible for the driving 
force behind dissolution. Absorption is internal, solid-phase, 
solvent accumulation, and may produce an initial swelling. 

Solid-phase activity competes with miscible displacement for 
leaching status. The leaching process is initially dominated by 
dissolution, followed by a progressive miscible drive. Time scaling 
dissolution with displacement is crucial for simulation. Simulation 
predicts increased leaching efficiency for solvents with higher 


adsorption coefficients. 


9.1. Single-Cell Dissolution 

The special case of uniform sphere radii in a liquid with initial 
concentration, Co; devoid of convection currents, constitutes the 
Single cell dissolution of bitumen. The fluid pressure is constant, and 
dispersion rapid, such that spatial dependence disappears leaving the 
State vector eee) with time dependence only. The system void of 


solvent sources can be modelled by the radial dissolution equation 


coupled with the adsorption equation. 
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ic 
r) : vA 
c* ae (roe) +f F(R (t-t)ar} - Hie NS ies Ci se = 0 
An S220 Z 
-—-Nr) — + 4nmr N_ [a(C-C* = 
(o, on Mp 7 ot p pL ( J Q Seth) 
By introducing the dimensionless variables, 
Dt 
ia een 
0 
3 
Thetis Beds62) 
0 
Fy ee 
4tr 
the transformed dissolution and adsorption equations have the 
respective dimensionless forms, 
2 3 “a: 
—C* _ (* 2 ee 
Tq eee 2 | C at, ee ry4(t)FA(ty tac} 0 
Dik 5) 
2, 9C Z Se be 
[1 - (1-S)) ra] at, + 3(1-S,) rgp, (C-C J,= 0 


where So is the initial liquid saturation and the Damkohler number, 


Nna , is the characteristic ratio of the solvent adsorption rate from 


the liquid phase to the solvent diffusion rate into the Piterior OL 


the solid. 
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The state vector ay determines the phase behaviour of the 
closed system for a given initial liquid saturation. Dissolution 
kinetics are governed by two mechanisms — solvent absorption and 
solvent adsorption. 

The accumulation term of the radial dissolution equation 
determines solvent absorption from the liquid phase with the critical 
concentration, C*, as the capacitance. In the absence of adsorption 
(Ny, = 9); absorption will yield an increase in solid mass as a result 
of internal solvent diffusion bounded by the critical SO sh ets) The 
swelling characteristics of solvent absorption are shown in Figure 5.1-1l, 
as a function of the capacitance. The saturation asymptotes reflect a 
solvent saturated solid to full capacity. The constant concentration 
Signifies the absence of bitumen desorption. A maximum swell 


capacitance is determined from the lower bound to the equilibrium 


saturation, = ae 
Ct =f EHC 5.1(5) 


above which, pore space is not available to accomodate full solid 
expansion and the absorption process must be truncated at the zero 
liquid saturation. 

Adsorption is the driving mechanism for the dissolution process. 
As solvent is adsorbed into the solid, the local concentration is forced 
above the critical. That portion above the critical undergoes a non- 
reversible phase change resulting in the dissolution of solid material. 


As dissolution proceeds, the liquid concentration drops indicating 
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bitumen desorption (solvent adsorption). An equilibrium concentration 
is reached when the solid phase is completely dissolved. Figure 5.1-2 
illustrates the dissolution curves as a function of the Damkohler 
number in the absence of absorption (C*=0). 

When the absorption rate exceeds the adsorption rate, an initial 
swelling is followed by rapid dissolution. Rapid dissolution is a 
consequence of the adsorption-perturbed absorption-elevated concentration 
Tevel within the solid. Figures 5.1-3 (a) and (b) illustrate this 
phenomenon for non zero capacitance, C* = .2 and C* = .4, respectively. 

The linear adsorption assumption poses difficulties about the zero 
radius. Dissolution must terminate at the zero radius, hence the 
integro-differential equation is valid only while solid material exists. 
After the dissolution of all solid material, a zero radius must prevail. 
The overshoots observed in these figures are a response of the curve 
fitting routine to the imposed dissolution discontinuity. 

The slightest perturbation of a solvent saturated solid will 
create an avalanching phase change. Figure 5.1-4 illustrates such 
dissolution avalanches, for a low initial saturation with high 
Capacitance. Adsorption decreases rapidly as the liquid concentration 
approaches the critical, leaving absorption to elevate internal 
concentration levels to the critical. Furthér adsorption ''tips the 
balance" resulting in a massive dissolution drive. 

Single-cell dissolution involves a closed system without 
additional solvent sources, which may be physically interpretted as the 
dissolution of tar sand in a laboratory beaker. Leaching, however, 


incorporates a multitude of interconnected cells with some local solvent 
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source being dissipated by both convective and dispersive mechanisms. 


5.2. Damkohler Effects 

Solid phase activity competes with miscible displacement for 
solvent. Initially constricted flow within the bitumen saturated sands 
allows rapid solvent absorption into the solid. After a short period 
of initial swelling, a dissolution driven non-reversible phase change 
increases liquid pore space. The upstream constrictions are gradually 
relieved, and freshly injected solvent accumulates upstream while 
fluid transport carries dissolved bitumen downstream. A dissolution 
wave trails the sluggish displacement, creating a high porosity zone 
with increased mobility. Upstream dissolution is driven by adsorption 
while absorption elevates downstream solid-phase concentration levels. 
Downstream solid activity responds to the fresh upstream solvent with 
accelerated rapid dissolution ultimately equalizing the saturation 
gradient. A burst of dissolved bitumen becomes apparent in the 
effluent, and dissolution terminates with the completion of total phase 
inversion. A progressive miscible drive follows to displace the 
remaining bitumen in solution. 

Pressures peak high in response to initially low injectivities. 
As dissolution proceeds, the pressures drop exponentially raising the 
injectivities by several orders of magnitude. 

Many factors influence the leaching process. Time scaling 


dissolution with displacement is crucial for resolution of the saturation 


history. The dimensionless time ratio, 
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relates the dissolution time scale to the pore volumes of solvent 
injected, giving an upper bound for the simulation time necessary to 
resolve saturation histories. 

A laboratory scale five-spot model whose characteristics are 
listed in Table 5.2-1 is chosen as a basis for simulation, giving a 
dimensionless time ratio (t4/t = .3) sufficient to resolve dissolution 
within the first pore volume of injection. 

An initial solvent invasion saturation is chosen as S, = .05, 
such that a maximum swell capacitance can be fixed accordingly to 
C* = .05. Simulations were executed at different adsorption rates, 
showing the Damkéhler effect on the leaching process. Figure 5.2-1 


graphs predicted production data for (a) ee oD) Ny a2) and 
a a 


(ce) Np = 5. All three runs illustrate a sharp increase in effluent 
a 


bitumen as a result of initial dissolution, followed by a gentle 
decline due to displacement. The first peak represents a stabilized 
system where the miscible drive balances bitumen desorption. As 
dissolution progresses, increased upstream porosity causes solvent 
to accumulate in the liquid phase. This accumulation drives down- 
stream dissolution to completion, producing a second tipple in the 
effluent bitumen. 

As the desorption rate increases with Damkohler number, rapid 
bitumen accumulation in the effluent causes a lower relative displace- 


ment activity. This is reflected by the higher peaks in the effluent 
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Figure 5.2-1 : Leaching - Damkohler Effects on Production 
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bitumen observed with increased Damkohler number. Simultaneously, 
solid depletion is accelerated and the miscible drive advanced. This 
is apparent in the steeper slopes of the bitumen decline. Higher 
cumulative bitumen to solvent ratios with larger recoveries suggest 
increased leaching efficiency. Economic efficiency also depends on the 
bitumen to solvent slope, such that a termination time may be predicted 
when production no longer supports leaching economics. 

A (6x6) diagonal grid was sufficiently large to damp any gross 
oscillation generated from numerical instability. Initial adsorption 
was significantly large such that TCP weights could not be obtained, 
and the second-order correct central difference scheme had to be 
adopted. The concentration surface was sufficiently smooth and 
diffusion sufficiently large to dissipate any convective oscillations. 
An initial indentation of the concentration front, similar to mild grid 
orientation, developed about the inlet but gradually disappeared as the 
miscible drive progressively dominated. 

The grid size was chosen small in order to reduce the enormous 
time of evaluating the series-solution of the Stephan problem, and the 
subsequent calculation of pointwise convolution integrals. Convergence 
by successive substitution was slow while dissolution prevailed. 
Problems were encountered about the zero radius where dissolution about 


singular points were forced to completion. 
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CHAPTER SIX 


CONCLUSION 


An overall discussion summarizes both theoretical model and the 
numerical simulator. Limitations are established, and leaching 
viability addressed. 

Achieved objectives are recognized with further proposals 


directed towards future research. 


6.1. Discussion 

The nonequilibrium leaching model predicts instantaneous porosity 
and saturation-dependent adsorption, in comparison to previous 
investigations for solution mining, limited to small porosity change 
with first-order reaction. An integro-differential absorption term keeps 
track of the solvent accumulated in the solid over the total time 
elapsed. Together with first-order adsorption, the instantaneous 
saturation follows from material balance. The solid phase spatial 
concentration distribution is implicit in the analytical solution of the 
nonstationary Stephan problem, and only the critical surface concentration 
at the spherical solid-liquid interface is needed a priori. 

Symmetry is essential for predicting similar dissolution fronts. 
Non-spherical particles experiencing radial dissolution will require a 
separate relationship for the dynamic geometry of the dissolving particle 
surface. Boundary layer effects may distort radial dissolution with 


accelerated upstream erosion at the leading particle edge. 
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The model confines each particle to its respective pore for all 
times of its dissolution history. However, fluid inertia may carry the 
particle intact, to some new location downstream. 

Local uniformity of the particle, in conjunction with a number 
density, relates the microscopic dissolution of individual particles to 
the bulk solvent transfer. For a large variance in regional particle 
size, a statistical extension to the integro-differential process must 
provide the macroscopic transfer rate. 

First-order adsorption assumes a linear concentration drop across 
the solid-liquid interface, with the ambient fluid concentration approxi- 
mated by the local macroscopic average. The model assumes a constant 
adsorption coefficient, and rapid interstitial mixing to equalize local 
concentrations. In the absence of absorption, linear first-order 
adsorption predicts asymptotic saturation histories. When combined with 
absorption, dissolution predicts finite particle-life times with a marked 
discontinuity at the zero radius. The critical surface concentration 
determines the particle capacity for solvent absorption, and extremely 
sensitive systems may result in highly accelerated dissolution from an 
avalanching phase change, terminating the particle. 

The linearity of adsorption creates problems in the presence of 
nonlinear absorption, with a singular equation about the zero radius. 

The linear model inherently assumes the ambient fluid concentration to be 
greater than the critical solid concentration for adsorption to proceed 
in the positive sense. Furthermore, accumulation within the transition 
zone is neglected. A more rigorous approach might model the individual 


pore-particle system, accounting for the microscopic variations in fluid 
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solvent to predict nonlinear adsorption. Fluid solvent would accumulate 
in a laminar boundary layer to yield a sharp spike in the normal 
concentration profile at the solid liquid contact surface. The non- 
reversible phase change would liberate heat locally about the dissolution 
site to be subsequently dissipated by the streaming ambient solution. 

An initial liquid continuity is essential for the model to function 
properly. The initial conditions investigated establish solvent to 
uniformly saturate the medium throughout. Such uniform sweeps do not 
accurately represent breakthrough conditions, and without better 
knowledge of the fingering mechanisms involved, the fully invaded medium 
becomes a basic mathematical device for initialization. 

The system responds to the uniformity of the initial conditions 
by adjusting to a stabilized state, where dissolution equalizes displace- 
ment at every point in the flow field. The stabilized concentration and 
saturation surfaces reflect a natural equilibrium state, characteristic 
of the flow geometry. As dissolution drives saturations higher, the 
increase in liquid pore space allows displacement to progressively 
dominate. Competition for solvent between absorption and displacement 
trigger rippling oscillations throughout the concentration surface. 
Bursts of bitumen appear in the effluent as the ripples migrate down- 
stream. 

The extent of rippling depends not only on the physical system, 
but also on the choice of initial conditions. If one could initiate 
leaching from the moment of solvent injection, fingering displacement 
mechanisms would determine the course of diffusing species, and the non- 


uniformity of the etched fluid paths would enhance pulsating of the 
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effluent bitumen. 

The leaching model has a closed form solution approximated with 
finite differences. The leaching simulator developed uses a five-point 
scheme for simplicity of implementation. Numerical limitations at 
high viscosity ratios become apparent when convection dominates. The 
partial differential equations become locally hyperbolic, where the 
finite-difference solution is unable to uniquely track the self-sharpening 
front. Mobility and diffusivity nonlinearities are approximated down- 
Stream (and back in time) to yield realistic miscible displacement on 
a diagonal grid. The shape of the predicted concentration surfaces are 
well represented, with channel formation at the extreme viscosity ratios 
associated with tar sand. 

Miscible similarity is controlled by the Peclet number, depending 
on the injection rate, grid partition, grain size and viscosity ratio 
in two dimensional displacements. The Damkohler number governs dissolution 
kinetics with the critical concentration as the absorptive capacity. Time 
scaling dissolution with displacement is crucial for the resolution of 
saturation histories. Dissolution is scaled by a characteristic 
dimensionless time dependent on both molecular diffusion rate in the solid 
and initial particle radius. A dimensionless time ratio relates the 
dissolution time scale to the pore volumes of injected fluid to estimate 
the simulation time necessary to resolve particle lifetimes. 

Simulation predicts increasing leaching efficiency for solvents 
with higher adsorption coefficients. Hot solvent injection is a viable 
procedure to increase adsorptivity, and enhance dissolution. However, 


a thermal model must be developed including a heat diffusivity equation 
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to determine temperature distributions. Hot solvents may invade tar 
sands more evenly, dissipating heat to effectively increase downstream 


mobilities. 


Ooo. rurther Study 

The thesis has presented an isothermal model for the in situ 
solvent leaching of tar sand. A two-dimensional, two-phase, two- 
component compositional simulator was developed to approximate the 
differential system. Analysis of simulated results identified parametric 
scaling effects providing insight for large-scale operations. 

Future work may investigate the thermal-compositional effects 
in a three-dimensional, multicomponent system with multi-fluid phases. 
A nonlinear adsorption model may be developed. High-accuracy numerical 


techniques may be implemented. 
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THE ANALYTICAL SOLUTION TO THE RADIAL SPHERICAL 


DIFFUSION EQUATION WITH A TIME DEPENDENT BOUNDARY 


APPENDIX A 


The method of Laplace transform is used to obtain the analytical 


solution to: 


satisfying the initial and boundary 
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follow from the change of variable: 
C= rC. A(10) 


C* = r C* A(11) 


Applying the Laplace transform, 
Vic@-tilm=el Ce cat = C(z,s) A(12) 
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to equation A(6), 
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interchanging the order of integration and differention in the first 


term, 


: | 
p. ae tongentane He eter sat 90 A(14) 
Ss 2 0 0 


i eee ere, ech 210 A(15) 


By virtue of the initial condition A(7), the middle term 
vanishes and the Laplace transformed equation becomes the ordinary 


differential equation, 
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The particular solution to equation A(21) satisfying the boundary 
eonditions “A(17)\ andgeA(l8) “as: 
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The solution to A(6) satisfying the initial and boundary 


conditions A(7), A(8) and A(9) is the inverse Laplace transform of 


A(26). 
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Define A, 
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Using the change of variable A(10) and A(11) gives the 
analytical solution to the radial diffusion equation A(1l) satisfying 


the initial and boundary conditions A(2), A(3), A(4) and ACS). 
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APPENDIX B 


DIRECTIONAL TRUNCATION CANCELLATION PROCEDURE 


FOR CONVECTION DOMINANT DISPLACEMENTS 


The truncation cancellation procedure, TCP, described herein 
was developed by the author for direct solution, paralleling the scheme 
developed by Laumbach*/ for the alternating direction implicit proce- 
Apes (ADIP). Unlike the ADIP version, the direct TCP involvesia 
transformation of the truncation error into velocity coordinates, where 
optimal weights suppress a dominant longitudinal error term. The 
analysis assumes convective dominance, where both adsorption and dif- 
fusion are relatively small. This approach deals with the differential 
form of the convective-diffusion equation and an extension to the 
integral form follows directly. 

The differential form of the convective-diffusion with adsorption 


(CDA) equation can be written in the form: 
\ ae i. 3C 2 vc. 
Ve (DevC) - V+ (UC) = ¢ 5 ae anh Flake ea | B(1) 


The CDA equation is discretized about mesh point, eee at the 
mid-time level, nts using central differences. Let A denote the 
central-difference operator, with (A, 24.) as the component spatial 
differences and A. as the temporal difference. Let ( ) denote the 


arithmetic time average relating the midtime values to the upper and 


lower time nodes. 
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The diffusion term is discretized to second-order accuracy in 


both space and time, 
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with the order of accuracy estimate, h. 
h = max (Ax,Ay,At) B(3) 


The adsorption term is time averaged to second-order accuracy. 
The component convective terms are discretized to second-order 


accuracy in both space and time, 
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to fourth-order. 
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The accumulation term is spatially weighted over the five-point 


difference star, 
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with E.. as the weighted truncation error of discretization also 


expressed to fourth-order. 
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The discretized CDA equation has the forn, 
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Ey and ES are the second-order errors of diffusion and adsorption, 


respectively, and the parameters T) and r, are the ratios of time 


step to component cell-dimensions. 
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As the rate of dissolution drops, the solid phase becomes 


virtually inert and miscible transport dominates. At high flow rates, 
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diffusion becomes small and convection dominates. In the limiting case, 
the CDA equation reduces to convection equation, 
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and the porosity becomes constant in time. With these assumptions, the 


third-order partials can be calculated directly from B(11). 
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Upon substitution of these third-order partials into B(9), the 


net discretization error becomes: 
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The weights Wy and Ww, are local to the immediate neighbourhood 


of the volume element about mesh point iG and are spatial averages 


with time dependence only. The differential operators are regrouped, 
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The error flux, ae is a second-order differential equation and 
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Evaluation of the dot products in B(16) gives the cartesian 


error components. 
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The error flux is dependent on both velocity magnitude and 
angle, such that a similarity transformation into velocity coordinates 
give rise to longitudinal and transverse components. In the velocity 
frame ane cear the velocity direction coincides with the principal 
longitudinal axis. The similarity transformation from cartesian to 


velocity frame involves a point rotation through a velocity angle, 6. 
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to give the error flux in velocity coordinates. 
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allows an isolation of the individual longitudinal and transverse error 


flux terms. 
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Termwise expansion of the component fluxes gives second-order 


differential expressions in the velocity frame. 
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The quantized flow within the volume element is assumed 


incompressible and irrotational, 
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such that the quantized velocity is constant about the immediate 
neighbourhood of the central mesh point. The quantized velocity is also 


assumed stationary about the mid-time node. 
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The partials are expanded, and then regrouped to yield: 
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The limiting convection equation in velocity coordinates, 


peau Se B(28) 


is differentiated with respect to the longitudinal coordinate, and with 
respect to time, to relate the second-order time derivative to the second- 


order spatial derivative. 
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Substitution into B(27) gives an expression for the quantized 


error flux, fully in terms of spatial partials. 
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Large concentration gradients may exist as a result of sharp 
solvent fronts propagating in the longitudinal direction. With 
transverse gradients small in the absence of diffusion, the dominant 
term of the second-order error flux is due to the longitudinal 
perturbation along the slopes of the propagating concentration front. 
The error flux, and consequently the divergence can be minimized by 
constraining the coefficients of the second-order longitudinal partials 


to vanish in both longitudinal and transverse components. 
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B(31) 


In this manner, the local errors from high convection rates are 
damped about the mesh point, and not allowed to grow with propagation. 
Both first-order numerical smearing and second-order longitudinal 
oscillations are minimized, but transverse perturbations may still 
force the second-order error to accumulate and its divergence effect 
the overall accuracy of the discretization. 


The system of constraints, B(31), is expanded in terms of 


component weights, 
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from an alternating direction implicit (ADI) procedure. 

The present analysis is valid when both diffusion and adsorption 
are negligable, and the optimal weights may differ considerably other- 
wise. As diffusion and adsorption become significant, their second- 
order errors, Ey and Ew effectively reduce the overall accuracy. A 
stability analysis determines the range of TCP applicability and gives 


the qualitative effects of larger diffusion and adsorption terms. 


STABILITY ANALYSIS 
The von Neumann method is applied to determine the stability 


criterion for the direct TCP discretization. The purely-convective 
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problem has the discrete form, 
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Cancelling the current central error from both sides simplifies 


the expression to: 
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In the absence of both velocity and porosity gradients, B(40) 


reduces to the complex equation, 


tim, = 0 B(41) 


with the real and imaginary parts, 
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determines this special case as weakly stable. 


The effects of a variable velocity field can be deduced from the 


a aaa | 
d ; ' a 
@ 4: ) Bee ; 
co 4 | v as 
aattifente xahte «sod mot? t9Ft9 tev 89 repairs oft gat 
#7, Ay. 
toe: aateeen — 2 
We . a. 
Lv ne is : 
‘ Le aS 
rp [-» of ? 
a : L 


barg YIieoror SHE v3 so at +, a. ons mt — 
ae 


161.8 ie 4 - oe 
(Pa 2 FTSEDRTo VI he 
: a * 


7 naw | suse YIRAGH sm bw os had " 


, ee By 7 
: i . : ' eer vite 
. : @ . 7 , } : : 


173 


first-order Taylor expansion about the mid-point velocity components in 


B(40). An additive real term, 


At Ue oy 
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implies adverse stability for a decreasing velocity field. 
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with phase dependent stability implications. 
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results in the addition of a real tern, 
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Re = Re, ies: A), + Reg] B(49) 


which effectively stabilizes the discretization. 
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Similarily, the inclusion of an adsorption term, 
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(Re, + Re +Re,+ Re.) + tm, + Im) 


A= : 
(Re, -Re -Re,-Re_) - t(im, - Im) 
0 Uu d a 0 ) B(53) 
Re, < 0 = STABLE 
The complex error equation of discretization will have the 
general form: 
z = Re + iIm = 0 B(54) 


To satisfy this equality, both real and imaginary parts must vanish 
simultaneously. For the special case of pure convection B(42), the 


system 
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to obtain the stability criterion for the weights. 
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This equality establishes the upper limit to the general case. 
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1 
as B(58) 


bounds the weights and fixes the maximum allowable time step for 


Stability. 


B(59) 


As Ay > ~, the discretization collapses to one-dimension, the 
velocity aligns with the x-direction and the y-weight vanishes. The 


nontrivial weight and the corresponding time step, 
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reduce to those cited by Laumbach?/ for the one-dimensional analysis. 
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APPENDIX C 


THE ANALYTICAL SOLUTION FOR THE ONE-DIMENSIONAL 


CONSTANT-COEFFICIENT CONVECTIVE-DIFFUSION EQUATION 


The method of Laplace transform is used to obtain the analytical 


solution to 


C = 0 rei, ‘eae ama 8 
C = Co —F=0, TS * 0 
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aes —E= 1, a ie | 


Applying the Laplace transforn, 


HiC(E, ct) = f eet ae = E(e,p) 
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Interchanging the order of integration and differentiation, 


: pd or C 
= ne fe Di dat es er Gdr|-"fP * pt o= dr = 0 
Z 0 T 


177 


C(1) 


C(2) 


C(3) 


C(4) 


C(5) 


C(6) 


C(7) 


: 
— 


4 : a 2 J , . 7 
“et % — a ay: a 


JAVOTaeaMia-ano “Sar OT wOITUIOG UOFTIIAMA om 
raed VOTREEAEA-RV Pra sero ast 


as 
Ce eS 


Se Cipne e i Sle wer arr A 
SF do oF bsen ra t? sosigsid to bolton sft : 
ane 
a {V's GES 
oe f.- ys ; a 


ey 7 = _ .< 


pé2 ThAGD atbaubd bite sabsint ort gniytetoae 


178 


and integrating the last term by parts, yields: 
~ pCO C(8) 
By virtue of the initial condition C(2), the third term of “C(8) 


vanishes and the Laplace transformed equation becomes the ordinary 


differential equation, 
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where A and B are constants to be determined by the boundary 
conditions C(3) and C(4). 


By defining the variable, 


Ce 4 a C(15) 


the general solution C(14) and its first-order spatial derivative can 


be written in the forms, 
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The particular solution to C(9) satisfying the boundary 


conditions: 6610) fand=-C (11) ies: 
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Co q+ e + 


ay sno 2 ie Crop Pea 
Pp JN. ”, / of/ 
i sys z ‘ped i oe as Nya’ 
ie me 
C22) 


The analytical solution to C(1) satisfying the initial and 
boundary conditions C(2), C(3), and C(4) is the inverse Laplace 


transform of C(22),.. 


Cte | CEs) C (23) 


By defining the variables, 


rages C (24) 
St*=sp + a2 =/G C(Z5) 


and applying the translation property of Laplace transform, 
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The first term of the series is recognized as the solution to 
the semi-infinite case and the physical significance of the successive 
terms constitute superimposed diffusion reflections at the outlet 
boundary (&=1). It can be shown that the series C(32) converges 


rapidly for Peclet numbers (N and times (t) of interest, and can 


Pe) 
be accurately approximated by the first three terms. 
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G(34))- into - C32), the analytical solution to the one-dimensional 
convective-diffusion equation with constant coefficients in terms of 


the Peclet number is: 
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APPENDIX D 


THE OPERATOR COMPACT IMPLICIT METHOD 


The operator compact implicit method>° seeks a linear relationship 
between the operator and the corresponding concentrations on the m-point 


difference star. For the one-dimensional constant coefficient case, 


3° 3 
EG) = pa yi 2 = a Je D(1) 


a; h(C,_1) + ayb(Cy) * ay 4y hCG 47) 
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So megane pe)” 0 
where the q's and r's_ are unknown weighting coefficients to be 
determined on the basis of the highest-order accuracy possible on the 
m-point difference-star. 

A linear combination of the concentration and their corresponding 
operators at the peripheral mesh points using fourth-order Taylor series 


results in the expression: 
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where the y-coefficients arise from the regrouping of Taylor terms. 
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forces the right-hand side of D(3) to be a linear combination of the 


concentration and operator at the central mesh point, 
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equation D(2) is obtained directly from D(6), and can be written in 


the compact form: 
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where S and I are the shift and identity operators, respectively. 


The compact operator can be isolated as 
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if the inverse of Q exists, and subsequently equated to the 


accumulation term, 


ac, 
ECW) ap D(11) 


Using the Crank-Nicolson temporal discretization gives the 
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The significance of To is observed when the boundary conditions 
are applied on the boundary-stars (44741 it 2).." For these difference- 
stars, the first-order spatial derivative is specified as a boundary 
condition and consequently the first-order y-coefficient need not be 
constrained. In order to obtain four equations in four unknowns, the 


Taylor series expansions are applied to fifth-order, such that 
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for inlet and outlet points, respectively. 
The invertibility of Q is discussed by Ciment, Leventhal and 
; 3 3 : é 5 eye e : ; 
Weinberg eh and for an interior point, the invertibility criterion 


depends on the cell-Peclet number, 


u Ax 


< 412 3.464 D(19) 
XX 
Extension to two dimensions for ADIP involves explicit 
evaluation of Onre however, for direct methods no explicit 
evaluation of mid-time operators are required. A complete two- 
dimensional operator including second-order mixed partials requires a 


nine-point (m=9) difference-star for direct solution methods. The 
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invertibility criterion limits the lower bound of the grid size necessary. 
Since 2(m-1) equations must be solved simultaneously at every point 
in the mesh-centered grid, the amount of work necessary to generate a 


solution becomes phenomenal in comparison to other methods. 
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